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Abstract 

The main theme of this dissertation is the study of the lattice points in a rational 
convex polyhedron and their encoding in terms of Barvinok's short rational functions. 
The first part of this thesis looks into theoretical applications of these rational func- 
tions to Optimization, Statistics, and Computational Algebra. The main theorem 
on Chapter |21 concerns the computation of the toric ideal I a of an integral n x d 
matrix A. We encode the binomials belonging to the toric ideal Ia associated with 
A using Barvinok's rational functions. If we fix d and n, this representation allows 
us to compute a universal Grobner basis and the reduced Grobner basis of the ideal 
Ia, with respect to any term order, in polynomial time. We also derive a polynomial 
time algorithm for normal form computations which replaces in this new encoding 
the usual reductions of the division algorithm. Chapter 01 presents three ways to 
use Barvinok's rational functions to solve Integer Programs: The {A,b,c)-test set 
algorithm, the Barvinok's binary search algorithm, and the digging algorithm. 

The second part of the thesis is experimental and consists mainly of the software 
package LattE, the first implementation of Barvinok's algorithm to compute short 
rational functions which encode the lattice points in a rational convex polytope. In 
Chapter 13 we report on experiments with families of well-known rational polytopes: 
multiway contingency tables, knapsack type problems, and rational polygons, and 
we present formulas for the Ehrhart quasi-polynomials of several hypersimplices and 
truncations of cubes. We also developed a new algorithm, the homogenized Barvinok's 
algorithm to compute the generating function for a rational polytope. We showed that 
it runs in polynomial time in fixed dimension. With the homogenized Barvinok's 
algorithm, we obtained new combinatorial formulas: the generating function for the 
number of 5 x 5 magic squares and the generating function for the number of 3x3x3x3 
magic cubes as rational functions. 



CHAPTER 1. Introduction 



Chapter 1 



Introduction 



1.1 Basic notations and definitions 

In this section, we will recall basic notations and definitions. Let {xi,X2, • • • ,Xm} be 
a finite set of points in M"'. A point 

m m 

X = 2^ OiiXi, where y, Q^j = 1 and aj > for z = 1, 2, . . . , m 

i=l i=l 

is called a convex combination of Xi, X2, . . . , Xm- Given two distinct points x, y G 
M*^, the set [x, y] = {ax + (1 — a)y : < a < 1} of all convex combinations of 
X and y is called the interval with endpoints x and y. A set C C M*^ is called 
convex, provided [x, y] G C for any two x, y G C. For C C M'^, the set of all convex 
combinations of points from C is called the convex hull of C and denoted conv(C). 
Let Ai, ^2, . . . , An e W^ and let 6i, 62, • • • , ^n e M. Then the set 

P := {x e R'^ : Ai ■ X < h ioT i = 1,2, ... ,n} 

is called a polyhedron. The convex hull of a finite set of points in M'^ is called 
a polytope and the ^yevl- Minkowski Theorem says that a polytope is a bounded 



polyhedron (jSchriiver 



19861 ) 
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A finite set of points {xi,X2, ■ ■ ■ ,Xk} C M*^ is afRnely independent if VAj G M, 

E •=! ^jXj = 0, E •=! A, = ^ A, = Vj = 1, 2, . . . , fc. 

A (d — 1) dimensional afiine set in Mf^ is called a hyperplane and every hyperplane 

can be represented a.s {x E M.'^ : ax = b, a E M.'^, a j^ 0, b E M}. a is called a normal 

vector of this hyperplane. 

Let iJ := {x G M'' : /i ■ X < /?}, where h E R'^, h ^ 0, and /3 G M, be an afRne 

half space. Then ii P C H and P n {x E R"^ : h ■ x = f3} ^ ^, then H is called 

a supporting hyperplane of P. A subset F of P is called a face if F = P or 

F = Pn H, where if is a supporting hyperplane. If a face F is minimal with respect 

to inclusion and F contains only a point, then F is called a vertex. 

Let V{P) be the set of all vertices of P. If any points in V{P) are in Z*^ then P is 

called an integral polyhedron. If any points in V{P) are in Q*^ then P is called a 

rational polyhedron. 

Now we will define our main tool, a generating function of a polyhedron. Let 

P C M"' be a polyhedron and let Z'^ C M^ be the integer lattice. For an integral point 

m = {nil, '"^2; • • • ! 'nT'd) £ '^'^j "we can write the monomial 

^ •— ^1 ^2 • • • ^d 

in d complex variables, zi, Z2, ■ ■ ■ , Zd- The generating function /(P, z) of a polyhedron 
P is the sum of monomials such that: 

fiP,z)= Y. Z-. (1.1) 

For example, consider the integral quadrilateral shown in Figure fTTD with the vertices 
Vi = (0,0), V2 = (5,0), V3 = (4,2), and 14 = (0,2). Then we have the generating 
function f{P,z) such that: 

/(P, Z) = Zi^ + 2:1^2 + Zi^ + 2:1^2:22 + 2:2^1^ + Zi^ + Zi^Z2^ + 2;22;i2 + Zi^ + Zi^Z2^ + 2;i2;2 + 

zi + 2;i2;2^ + 2:2^ + Z2 + I. 
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Figure 1.1: An example for the generating function 



One notices that the multivariate generating function f{P, z) has exponentially many 
monomials even though we fixed the dimension. So one might ask if it is possible 
to encode /(P, z) in a "short" way. In 1994, A. Barvinok showed an algorithm 
that counts the lattice points inside P in polynomial time when rf is a constant 



( Barvinok . 



1994^ . The input for this algorithm is the binary encoding of the integers 



Aij and 6j, and the output is a short formula for the multivariate generating function 
/(P, z) = XlaePnZ'' '^"- This long polynomial with exponentially many monomials is 
encoded as a short sum of rational functions in the form 



fiP, ^) 



E 



'1.21 



ie/ 



;i -^^1.^(1 - z""^'^) . . . {I - z'^d'^ 

where Uj, Ci^j, C2,i, . . . , Cd^i € TJ^ and where / is a polynomial sized index set. 
We call this short sum of rational functions of the form p.2|) Barvinok's rational 
function for the generating function f{P,z). For brevity, we also call it a short 
rational function for the generating function f{P,z). For example, suppose we 
have the polytope in Figure fTTTl Then we can write: 

/(P, Z) = Zi^ + 2:1^2:2 + Zi^ + Zi^Z2^ + Z2Zi^ + Zi^ + 21^22^ + 22^1^ + ^1^ + ^1^^2^ + ^1^2 + 

zi + 2:1 2:2^ + 2:2^ + ^2 + 1 



+ 



Z^Z^ 

[X-z-'zl-^^X-z-'Y 
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Here is another example to clarify Barvinok's rational functions. Suppose we have a 
tetrahedron P with vertices v^ = (0,0,0), V2 = (1000000,0,0), V3 = (0,1000000,0), 
and ^4 = (0, 0, 1000000) in Figure 11.21 Then we have the multivariate generating 
function f{P, z) which has 166, 667, 666, 668, 500, 001 monomials. However, if we use 
Barvinok's rational functions, we can represent all these monomials using a "small" 
encoding: 



f{P.^) 



+ 



^1000000 

^1 



1 - zi)(i - Z2){i - 2:3) (1 - ^r'^2)(i - ^r'^3)(i - ^r'; 



4000000 



4000000 



(1 - Z,Z:,'){1 - Z:,'Z,){1 - Z^') (1 - Z,Z^')il - Z2Z^')il - z^'Y 




Figure 1.2: A tetrahedron example for the generating function 



One might ask how one can compute Barvinok's rational functions for the input 
polyhedron. Th e following theorem tells us that there is an algorithm created by 
Barvinokl (|199J) to compute Barvinok's rational functions from the input polyhedron 
in polynomial time in fixed dimension. We will describe a variation of the algorithm 
in Chapter 0J 
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Theorem 1.1. Warvinoli. \l994- Theorem 5.4) Fix the dimension d. Then there 
exists a polynomial time algorithm which for a given rational polyhedron P C M'', 
computes f{P, z) in the form of M.^) in polynomial time. 

We do not want to expand Barvinok's rational functions because expanding them 
causes exponential complexity. So, if we want to perform operations on sets via gen- 
erating functions, such as taking unions, intersections, projections, and complements, 
we want to do it directly with Barvinok's rational functions without expanding them. 
The Hadamard product of Laurent power series is a very useful tool for Boolean 
operations on sets via Barvinok's rational functions. 

Definition 1.2. Let gi and g2 be Laurent power series in z & C^ such that gi{z) = 
Ylim&zd ^mZ'^ (ind (72 (-z) = J2mezd ^mZ^ ■ Then the Hadamard product g = gi * g2 is 
the power series such that: 

9{z) = ^ ambmz"^. 

Hadamard products of Laurent power series are one of the most important tools to 
prove theorems in this thesis. They are used for taking unions of sets, intersections of 
sets, and set difference via short rational functions without expanding them. We will 
show how to compute the Hadamard product of Laurent power series gi and (72 given 
in the form of rational functions. Let pi,P2, 0,11, ■ ■ ■ , dik £ ^'^ and 021, . . . , a2fc G ^'^• 
Suppose we are given the Laurent power series gi and (72 in the form: 

9i = -, ^ ; r and 02 = 7 ; ; r. (1.3) 

^ (l-2'^ii)...(l -^'^ife) ^ {1 - z''^^) . . . {1 - z^^k.) ^ ^ 

Here is an outline of the algorithm to take the Hadamard product of two Laurent 

power series via Barvinok's rational functions. 



Algorithm 1.3. Warvinok and Wood; 



mm, Lemma 3.4) 



Input.- Laurent power series gi and (72 iiT- the form of l\l.S\) . 
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Output; The Hadamard product gi*g2 of gi and g2 in the form of a rational function. 
Step l; If aij > or a2j > 0, then apply the identity: 



to reverse the direction of aij or a2j. 

Step 2: Let P C M^^ = {(^^^ 6, • • • , 6fc) : {i G M /or z = 1, 2, . . . , 2A;} &e a rational 

polyhedron defined by the equation: 

Pi + 6«ii H \- ikaik = P2 + ^k+ia2i H h 6fca2fc, 

anc? inequalities 

^i>0 for I = 1,2,..., 2k. 

Step 3; Using Theorem M . IV we compute: 



fip,x) = j2±- 



for Ui,Vij G Z^'^. 

Step 5; Apply the monomial substitution $ : C^*"' -^ C^ such that: 

to the function f{P,x). 

Step 4; Return gi * g2 = z^^$(/(P, x)). 

1.2 Computer Algebra and applications to Statistics 

One focus of this thesis is applying Barvinok's rational functions to Statistics and 
Mixed Integer Programming. First we consider the connection between Computa- 
tional Algebra and contingency tables in Statistics. 
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(1.4) 



Definition 1.4. A s-table of size (ni, . . . ,ns) is an array of non-negative integers 
"^ = ('yii....,ii,); 1 ^ "^i < iT-j- For < L < s, an L-marginal of v is any of the Q) 
possible L-tables obtained by summing the entries over all but L indices. 

Example 1.5. Consider a 3-table X = (xijk) of size {m, n, p), where m, n, andp are 
natural numbers. Let the integral matrices Mi = (ajk), M2 = (bik), and M3 = (qj) 
be 2-marginals of X, where Mi, M2, and M3 are integral matrices of type n x p, 
mxp, and mxn respectively. Then, a 3-table X = (xijk) of size {m, n, p) with given 
marginals satisfies the system of equations and inequalities: 

YJLi^ijk = ajk, U = 1,2, ...,n, k= 1,2, ...,p), 
Y.'^^iXijk = bik, {i = l,2,...,m, k = l,2,...,p), 

2_^k=l -Xijk Cjj; Y' i,Z, ...,?Tl, ] i,Z,...,?T.J, 

Xijk > 0, (z = l,2,...,m, j = l,2,...,n, k = 1,2, ...,p). 
Such tables appear naturally in Statistics and Operations Research under various 
names such as multi-way contingency tables, or tabular data. We consider the table 
counting problem and table sampling problem: 

Problem 1.6. (Table counting problem) 

Given a prescribed collection of marginals, how many d-tables are there that share 

these marginals? 

Problem 1.7. (Table sampling problem) 

Given a prescribed collection of marginals, generate typical tables that share these 

marginals. 

The table counting problem and table sampling problem have several applications in 
statistical anal ysis, in particular for independence testing, and have been the focus of 



1998; 



much research (lAnderson and Fienberd. 



Dobra and SuUivant 



2002; 



2001' 



RapalL 



I 



De Loera and Onn 



2002; 



Diaconis and Sturmfels . 



2003|). Given a specified collection of 



marginals for d-tables of size (rii, . . . , n^) (possibly together with specified lower and 
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upper bounds on some of the table entries) the associated multi-index transportation 
polytope is the set of all non-negative real valued arrays satisfying the given marginals 
and entry bounds specified in the system of equations and inequalities, such as for- 
mulas ()1.4|1 for a 3-table given in Example 11.51 The counting problem is the same as 
counting the number of integer points in the associated multi-index transportation 
polytope. 

In this thesis, one of the main tools to solve the table counting problems and table 
sampling problems is Computational Algebra. We consider a special ideal in the 
multivariate polynomial ring, namely a toric ideal I a associate to the given integral 
matrix A. We compute the Grobner basis associate to the toric ideal I a. Then we 
apply Grobner bases of the toric ideal I a to the table counting problem and the table 
sampling problem . H ere, we w o uld li ke to remind the reader of some definitions. 



Cox et al 



( 19971 ) and ISturmfeld ()1996[ ) are very good references for details. Let us 
denote Z^ := {x G Z'^ : x > 0} and Z+ := {x G Z : x > 0}. 

Definition 1.8. Let K be any field and let K[x\ = K[xi, X2, . . . , x^] be the polynomial 
ring in d indeterminates. A monomial is a product of powers of variables in K[x\, 
i.e. x"^X2^ . . ■x'^'^ , where ai, 0^2, . . . , a^ G Z4.. 

Definition 1.9. Let K be any field and let K[x\ = K[xi, X2, . . . , x^] be the polynomial 
ring in d indeterminates. Let I C K[x\. Then we call I an ideal if it satisfies the 
following: 

•f + geIforallf,geI. 

• af & I for all f E I and all a G K[x\. 



Note that by Hilbert basis theorem (jCox et al. 



19971 Chapter 2, section 5, Theorem 



4) every ideal in K[x\ is generated by finitely many elements in K[x\. 

Definition 1.10. Let ^ be a total order on Z^. We call -< a term order if it 

satisfies the following: 
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• For any a, j3, S E Z'^, a-<P-^a + S-<(3 + 6. 

• For any a G Z'^\{0}, ^ a. 

A term order on Z^ gives a term order on the monomials of K[x] by setting a bijection 
map from Z^ to S* C K[x], where S is the set of all monomials in K[x], such that 
(ai, a2, ..., ad) -* x^^x^^ ■■■x'^d- 

Definition 1.11. (The lexicographic term ordering) Let a, (3 E Z^. We say a -<iex P 
if the left most non-zero entry of a — (3 E 7/ is negative. We write x" -<iex x^ if 

a ~<lex P- 

For example, if we have (3, 2, 7) and (3, 5, 2) in M^, then we have (3, 2, 7) — (3, 5, 2) = 
(0,-3,5). So, (3,2,7) -<iex (3,5,2) and xfx^xl -<iex xfxlxl. One notices that the 
lexicographic term ordering is a term order. 

We can also define a term order from a vector c, as we described in Definition 11.101 
by the following method: we make this vector c into a term order -<c such that for 
all a, f3 e Z^, a ^e /5 if 

• €■ a < c ■ P or 

• ca = c/3 and a -<iex P- 

For example, suppose c = (1, 0, 2) and if we have (3, 2, 7) and (3, 5, 2) in M^, then we 
have (1, 0, 2) • (3, 2, 7) = 17 and (1, 0, 2) • (3, 5, 2) = 13. So, since (1, 0, 2) ■ (3, 5, 2) < 
(1, 0, 2) ■ (3, 2, 7), we have (3, 5, 2) ^^ (3, 2, 7) and xfx^xj -<c x\xlxl. 
In general, any term order is defined by a d x d integral matrix W . We represent 
a term order -< on monom ials in xi, . . . ,Xrf by an integral d x (i-matrix W as in 



fJMora and Robbiano . 



19981). Two monomials satisfy x° -< x^ if and only if Wa is 
lexicographically smaller than Wp. In other words, if wi, . . . ,Wd denote the rows of 
W, there is some j G {1, . . . ,d} such that Wia = WiP for i < j, and Wja < WjP. For 
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example, W = Id describes the lexicographic term ordering. We will denote by -<w 
the term order defined by W. 

Definition 1.12. Let K he any field and let K[x\ = K[xi,X2, ■ ■ ■ ,Xd] be the poly- 
nomial ring in d indeterminates . Given a term order -<, every non-zero polynomial 
f G K[x] has a unique initial monomial, denoted in^{f). If I is an ideal in K[x], 
then its initial ideal is the monomial ideal 

in^{I) ■.=<in^{f) : f el> . 

The monomials which do not lie in in^{I) are called standard monomials. A finite 
subset G d I is called a Grobner basis for I with respect to -< ifin^{I) is generated 
by {in^{g) : g G G}. A Grobner basis is called reduced if for any two distinct 
elements g,g & G, no terms of g is divisible by in^{g). 



Proposition 1.13. XCox et all \199'^. Proposition 1, Chapter 6) 

Let G := {(71, 5'2, • • • , fi'fc} be a Grobner basis for an ideal I C K[x\ and let f G /('[x]. 

Then there exists a unique r & K[x] such that: 

• No term of r is divisible by any of leading term of gi, for all i = 1,2, ... ,k. 

• There is g & I such that f = g + r. 

In particular r is the remainder on division of f by G, and r is unique no matter how 
the elements of G are listed when using the division algorithm. 

The remainder r for / G K[x] is called the normal form of /. Note that the re- 
duced Grobner basis is unique. This thesis concentrates in a special kind of ideals / 
in K[x] = K[xi,X2, . . . ,Xd\, which are called toric ideals. Toric ideals find applica- 
tions in Integer P rogramming, Computational Algebra, and Computational Statistics 



( Sturmfels 



19961). 
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Definition 1.14. Fix a subset A = {ai,a2, ■ ■ ■ ,ad} o/Z". 

Each vector ai is identified with a monomial in the Laurent polynomial ring K[±t] :- 



K[t,t^... 
mial map 



TA-\t 



1 +-2 



t "^J . Consider the homomorphism induced by the mono- 



vr: K[x] -^ K[±t], Xi -^ t°\ 
Then the kernel of the homomorphism -if is called the toric ideal of A. 

The following lemma describes the set of generators of a toric ideal I a associated to 
the integral matrix A. 



Lemma 1.15. XSturmfela . \l99d. Lemma 4-i) The toric ideal I a is spanned as a 
K-vector space by the set of binomials 

{x" -x" : u,ve Z^, Au = Av}. 

The main theorem on Chapter |21 concerns a new way to compute the toric ideal I a 
of the integral matrix A. 

Now we are ready to discuss applications of Computational Algebra to Computational 
Statistics. As we mentioned earlier, a toric ideal Ia and the Grobner basis associated 



19961 ). Here we would 



to I A find applications to Computational Statistics (jSturmfels . 

like to discuss how we can apply Grobner bases to solve the table counting and table 

sampling problems. First of all, we will remind the reader of the definit i on of Markov 



bases associate to the given integral matrix A (JDiaconis and Sturmfel 



19981) . 



Definition 1.16. Let P = {x e R'^ : Ax = b,x > 0} j^ 0, where A e Z"^^ and 
b G Z", and let M be a finite set such that M C {a; G Z"' : Ax = 0}. Then we define 
the graph Gb such that: 

• Nodes of Gb are lattice points inside P. 



Draw a undirected edge between a node u and a node v if and only if u — v G M . 
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Then we call M a Markov basis of the toric ideal associate to a matrix A if G^ is 
connected for all b with P ^ ^. If M is minimal with respect to inclusion, then we 
call M a minimal Markov basis. 

Note that, in general, a minimal Markov basis is not necessarily unique. A Markov 



basis c an be used for randomlv sampling data and random wa^ 



tables (JDiaconis and GaneoUi . 



1995 



Diaconis and Sturmfels 



ks on contingency 



19981). We will describe 



the Monte Carlo Markov Chain algorithm which uses Markov bases to create random 
walks on contingency tables. We can also define a Grobner basis using a graph Gb- 



Lemma 1.17. kSturmfelA . \l99f\. Theorem 5.5) Let P = {x ^W^ : Ax = b,x > Q} ^^ 

0, where A e Z"""^ and fe e Z". Let M be a finite set such that M (Z {x eZ'^ : Ax = Q} 
and let -< be any term order on 'H'^ . Then we define the graph Gb such that: 

• Nodes of Gb are lattice points inside P. 

• Draw a directed edge between a node u and a node v if and only if u -< v for 

u-veM. 

If Gb is acyclic and has a unique sink for all b with P ^ 0, then M is a Crobner basis 
for a toric ideal associate to a matrix A with respect to -< . 

Notice that if M is a Grobner basis then this implies M is a Markov basis because 
if we have an acyclic directed graph with a unique sink, then it has to be connected. 
However, note that not all Markov bases are Grobner bases. 

Remark: Grobner bases provide a way to generate Markov bases for a wide variety 
of problems where no natural set of moves were known. 

Example 1.18. Suppose we have 2x3 tables with given marginals. 
There are 19 tables with these marginals for 2x3 tables in Table Wl\ 
Up to signs, there are 3 elements in the Markov basis. 
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? ? ? 


? ? ? 


? ? ? 


6 










? ? ? 


? ? ? 


? ? ? 


6 








Total 


4 
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4 





Table 1.1: 2x3 tables with 1-marginals. 
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Figure 1.3: Markov basis elements for 2 x 3 tables 



Figure [7"^ gives a connected graph for a Markov basis for 2x3 tables. An element of 
the Markov basis is a undirected edge between integral points in the polytope. 

From here on, we focus on Markov bases for contingency tables. Why do we care about 
Markov bases for contingency tables? We care about them because using Markov 
bases, we c an estimate the number o f table s b y Monte Carlo Markov Cha i n (MC MC) 



algorithm (JDiaconis and Sturmfelsl (|l998l )). iDiaconis and Saloff-Costd ()1995[ ) also 
showed that the rate of convergence for MCMC is 6"^, where 6 is the diameter of 
the graph Gb- The outline of MCMC algorithm is the following: 



Algorithm 1.19. iDiaconis and Sturmfek . 



19981. Lemma 2. 1) 



(Random walk on a graph) 

Input: A Markov basis M of a graph Gb defined by the set of contingency tables with 

given marginals and an initial node /o in Gb. 

Output: A sample from the hypergeometric distribution a{-) on Gb. 
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Figure 1.4: A connected graph Gb f^or given h. 
1. Set f := /o and sei t = 0. 
^. While ft < S^) 

• Choose u & M uniformly and a sign e = ±1 with probability 1/2 each 
independently from u. 

• If f+eu > then move the chain from f to f + eu with probability Tain{a{f + 
eu)/a{f), 1}. // not, stay at f . Set t = t + 1. 



Diaconis and Sturmfelsl ()l998( ) showed that this random walk on the graph Gb is a 



connected, reversible, aperiodic Markov chain on Gh which converges to the hyperge- 



ometric distribution (^Diaconis and Sturmfeb 



19981 . Lemma 2.1). With the random 



walk on the graph Gb, we apply it to approximating the number of lattice points in a 
convex rational polytope P. The idea for approximating the number of lattice points 
in a convex polytope P is the following: 

Suppose we take a sequence of draws H := {ni, ?i2 • • • , Um} randomly from the uniform 
distribution over P. Let p{ni) = 1/\P\ be the uniform distribution over P. If we can 
simulate a lattice point n^ & P from a distribution g(-), where q(t) > for all t & P, 
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then we have 



E\ 



h=Y.^A^) = \p\- 



lit) 



teP 



lit) 



Hence by the Strong Law of Large Number (JDurrett . 
estimation of I PI is: 



2nnd (7.1) on page 56), the 



1 "" 1 
m j^ q{ni) 



1.3 Algorithms for counting 



Enumerating the lattice points in a given polytope and counting the number of lattice 



points in a given polytope is very useful to Computational Statistics (IDe Loera and Onn 



2002 



2002 



Diaconis 



Rapallo 



and 



GangoUi . 



1995; 



Diaconis and Sturmfels 



1998 : 



Dobra and SuUivant . 



2003^ . So, in the next section, we will briefly discuss methods to count 
exactly the number of lattice points in a given convex polyhedron. Barvinok's method 
will be fully discussed in Chapter HI 

Using the multivariate generating function f{P,z) for a polytope P, we can count 
the number of lattice points inside P. In fact, the number of lattice points inside P 
is/(P,(l,l,...,l)). 

Example: Let P be the quadrangle with vertices Vi = (0, 0), V2 = (5, 0), V3 = (4, 2), 
and 1/4 = (0,2). 
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Then we have: 
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/(P, Z) = Zi^ + Zi^Z2 + Zi^ + Zi^Z2^ + 2:2-21^ + Zi^ + Zi^Z2^ + 2:2-21^ + ^1^ + Zi^ Z2'^ + 2^1 2^2 + 
-21 + -21-22^ + Z2^ + 2:2 + 1. 

If we substitute zi = 1 and Z2 = 1 iti f{P, z), then we have f{P, (1, 1)) = 16, which is 

also the number of lattice points inside P. 

We can use Barvinok's rational functions to count the number of lattice points inside a 

polytope. Notice that, unfortunately, the point (2:1, 2:2, . . . , 2;^) = (1, 1, . . . , 1) is a pole 

of Barvinok's rational functions. Thus, we cannot directly substitute (2:1, 2:2, . . . , 2;^;) = 

(1,1,..., 1) into /(P, z). Instead, we compute lim2^(i ... 1) /(P, z), which is the number 

of lattice points inside P. In this thesis, we apply the residue calculus to compute 

lim^^(i 1) /(P, z). We will discuss how to compute this limit via the residue calculus 

in Section mm 

Severa l "ana l ytic" a lgorithms have been proposed by many authors 



(2002) 



19601 ) 



Beckl(l2003[): 



Lasserre and ZeronI (J2002i ) 



Pemantle and Wilson 



mm. 



Baldoni-Silva and Vergne 



Lasserre and Zeron 



20031 ) 



MacMahon 



A couple of these methods have been imple- 
mented and appear as the fastest for unimodular polyhedra. However, only Barvi- 
nok's method has been implemented for arbitrary rational polytopes. Consider, for 
example. Beck's method: let Mi be the columns of the matrix M. We can interpret 
P{M, b)n7j'^ as the Taylor coefficient of z^ for the function H^^^- — to--. One approach 
to obta in the partic ular coefficient is to use the residue theorem. For example, it was 



seen m 



Becd pOOOl ) that 



P{M, h) n If 



1 



(27ri) 



2 



-61-1 



z„ 



m — ^m 



(1 



z 



Ml' 



[l - Z^-' 



dz . 



ziNei J\zn 

Here < ei,...,em < 1 are different numbers such that we can expand all the 

rr- into the power series about 0. It is possible to do a partial fraction decom- 

1 — 2;^''^fe 

position of the integrand into a sum of simple fractions. This was done very s uccess- 



fully to carry out very hard computations regarding the Birkhoff polytopes (JBeck 



2003^ . Vergne and collaborators have recently developed a powerful general the- 
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ory about the multivariate ratio nal functions 11 



2002; 



Szenes and Vergne 



"i=i(i_^".) 



( Baldoni-Silva and Vergne 



2002) • Experimental results show that it is a very fast 



met hod for unimodular po 



son (Pemantle and Wilson 



ytop es (JBaldoni-Silva et al. . 



20031) ■ Pemantle and Wil- 



20031 ) have pursued an even more general computational 



theory of rational generating functions where the denominators are not necessarily 
products of linear forms. 



Recently, 



Lasserre and Zeron 



( 20031 ) introduced another method to enumerate the 
lattice points in a rational convex polyhedron. Suppose we have a rational convex 
polyhedron Py = {x E W^ : Ax = y, x > 0}, where A = (Aij) e Z"""^ and b G Z". 



Then we define the function 



/w- E 



1.5) 



xeP. 



where c G Z"^ is small enough so that f{z) is well defined. A tool which 
(120031 ) use is a generating function F : C" ^ C: 



Lasserre and Zeron 



z^F{z):=Y,fiy)^'- (1-6) 

yeZ" 

Note that the generating functions in (jl.Sp and (|1.6|) are different from Barvinok's. 
Let us define P* := {b e R"^ : b ■ x > 0, x e Pq} and L := {c G M'^ : -c > 
b, for some b G Pq}- Suppose Ai is the ith column of a matrix A. Then we have the 
following lemma. 



Lemma 1.20. \Lasserre and Zeron 



200: 



, Proposition 2.4) 



Let f and F be functions defined in M.5}) and il.6]) . and let c eT. Then: 



^^ = 1177 



1 



— e'-Zi 



^71 , 



fc = l V- - -1 

on the domain (|2;i|, . . . , |^„|) G {y G M" : y > 0, e'^'^x^'' <l,k = l,...,d}. 
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Definition 1.21. \Lasserre and Zeron I . \2003l. Definition 2.1) Let p G N satisfy n < 
p < d and let v = {ui, . . . , Up} C N be an ordered set with cardinality |z/| = p and 
I < Ui < ■ ■ ■ < i^p < d. Then 

1. V is said to he a basis of order p if the n x p suhmatrix Ay = [Aj,J . . . lA^, ] has 
the maximum rank. 

2. For n < p < d, let Jp := {u C {1, . . . , d}]!^ is a basis of order p}. 



Then 



Lasserre and Zeron 



( 20031 ) show how to invert the generating function F{z) 
in order to obtain the exact value of f{y). First, they determine an appropriate 
expansion of the generating function in the form: 



F(z) 



E 



Qa 



;i-7) 



where the coefficient Qo- : C" — i> C are rational functions with a finite Laurent series 



f3eZ",\\/3\\<M 

for a strictly positive integer M. Then they apply the following theorem to obtain 



Theorem 1.22. ^Lasserre and Zeron 



200i 



, Theorem 2.6) 



Let A G Z"^'^ he of maximal rank, let f he as in M.5]) with c G F. Assume that the 
generating function F in |7~^) satisfies \1.1\) and \1.^) . Then, 

o-eJn /3gZ",||/3||<M 



with 

E^y-P) 

where c^ = (c^i,. . •,c^„)- 



e^-^ if X := A-\y - f3) eW, 
otherwise; 
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The main point is that as soon as we have f{y) with sufficiently small c vector, if we 
send Cj — s> for i = 1,2, .. .d, then we can obtain the number of lattice points inside 
a convex rational polytope Py (if Py is a unbounded polyhedron, then this limit does 
not converge). 

1.4 Applications to Mixed Integer Programming 

Now we explain the connections between Barvinok's rational functions and integer 
programming problems. We consider an integer linear programming problem: 

Question 1.23. (Integer Programming) 

Suppose A G Z"^°', c G Z'^, and b G Z"'. We assume that the rank of A is n. Given 

a polyhedron P = {x ^W'' : Ax = h, x >Q}, we want to solve the following problem: 



(IP) maximize c ■ x subject to x G P, x G Z . 



These problems are ca 



lem is NP-hard (IKarp . 



l ed int eger progra mming problems and we know that this prob 



19721) . However, 



Lenstra 



( 1983| ) showed that if we fixed the 



dimension, we can solve (IP) in polynomial time. Originally, Barvinok's counting 
algorithm relied on H. Lenstra's polynomial tim e algorithm for Integer Programming 



in a fixed number of varia 
through 



es (ILenstra 



I 



1983^ . but shortly after Barvinok's break- 



Dver and KannanI ()1993l ) showed that this step can be replaced by a short- 



vector computation using the LLL algorithm. Therefore, using binary search, one can 
turn Barvinok's counting oracle into an algorithm that solves integer programming 
problems with a fixed number of variables in polynomial time (i.e. by counting the 
number of lattice points in P that satisfy c • a; > a, we can narrow the range for the 
maximum value of c ■ x, then we iteratively look for t he largest a where the count 
is no n-zero). This idea was proposed by Barvinok in ^Barvinok and Pommersheiml . 
19991 ). We call this IP algorithm the BBS algorithm. 
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To solve more general problems one might want to consider some variables as real 
numbers, instead of all integers. Thus one can set the following problems: 

Question 1.24. (Mixed Integer Programming) 

Suppose A G Z"^*^, c G Z"^, and b G Z". We assume that the rank of A is n. Given 

a polyhedron P = {x eM!^ : Ax = &, x > 0}, we want to solve the following problem: 

(MIP) maximize c ■ x subject to x ^ P, Xi & Tj if the index i G J G [d]. 

One notices that linear programming problems form a subset of (MIP) problems, i.e. 
J = % and also integer programming problems form a subset of (MIP) problems, i.e. 
J = [d]. 

As we mentioned earlier, we can define the term order from a cost vector c, as we 
described in Definition II. 101 From this term order -<c, as soon as we have the reduced 
Grobner basis with the term order ^c, we can solve the integer programming problem 
min{c ■ X : X E P \~\ U^}. A sketch of an algorithm is the following: 



Algorithm 1.25. \Sturmfela. \l99d. Algorithm 5.6) 

Input: A cost vector c G TJ^, a matrix A G Z"^*^, a vector b G 7/^ and a feasible 

solution vq E P f] Tf" , where P := {x E^^ : Ax = b, x > 0}. 

Output: An optimal solution and the optimal value of minimize c ■ x subject to 

xE PnZ'^. 

Step 1: Compute the Grobner basis with the term order -<c. 

Step 2: Compute the normal form x" of x'"° and return u and cu, which are an 

optimal solution and the optimal value, respectively. 

One notices that Algorithm II .251 outputs the optimal value and an optimal solution 
for minimization. Trivially if one wants to have the optimal value and an optimal 
solution for maximization such as (IP), one can set c = — c and apply Algorithm ll.251 
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Definition 1.26. Suppose we have a rational convex polyhedron P C M'^ and suppose 
we have an integer programming problem such that maximize c-x subject to x E P{~\l/. 
Also suppose a feasible solution Xq E P P^TJ^ is given. Then we call an integral vector 
t E Ij'^ an augmenting vector if c{xo + t) > cxq. A finite set which contains all 
augmenting vectors is called a test set. 

There has been considerable activity in the area of te st sets and au g menta ti on meth- 



ods (e .g. Graver, integral basis method, etc. See lAardal et al.l (|2002cr ): 



Thomas 



(120011 1). One notices from Algorithm 11.251 that the Grobner basis of a toric ideal I a 
associated to the integral matrix A with respect to a term order -<c is a test set for 
an integer programming problem min{c ■ x : Ax = b, x > 0, ,x E 1^'^}. 

In 2003, Lasserre observed a new method for solving integer programming problems 
using Barvinok's short rational functions, which is different from the BBS algorithm. 
We consider the integer programming problem maximize{c ■ x : Ax < b,x > 0,x E 
Z*^}, where c E Z,'^, A E Z'"^'^, and b E Z"\ This problem is equivalent to Problem 
11.231 because using Hermite normal form, we can project the polytope P down into 
a lower dimension until the dimension of P equals to the dimension of ambience 
space. We assume that the input system of inequalities Ax < b,x > defines a 
bounded polytope P C W^, such that P n Z*^ is nonempty. As before, all integer 
points are encoded as a short rational function /(P, z) in Equation (|1.2p for P, where 
the rational function is given in Barvinok's form. Remember that if we were to 
expand Equation (|1.2p into monomials (generally a very bad idea!) we would get 
f{P,z) = XlaGPnZ'* -^"- -^^^ ^ given c E Z*^, we make the substitution Zi = t^% 
Equation (|1.2j) yields a univariate rational function in t: 



f(p, t) = y ±—r^ . (1.9) 

The key observation is that if we make that substitution directly into the monomial 
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expansion of f{P, z), we have z°' —>■ f^". Moreover we would obtain the relation 

f{P,t)= Y^ r" = H*^ + (lower degree terms), (1.10) 

aePnZ'i 

where M is the optimal value of our Integer Program and where k counts the number 
of optimal integer solutions. Unfortunately, in practice, M and the number of lattice 
points in P may be huge and we need to avoid the monomial expansion step altogether. 
AH computatio ns have to be done by manipulating short rational functions in ()1.9|) . 
Lasserrd (J200J) suggested the following approach: for i G /, define sets rii by r/j = 
{j G {1, ...,d} : c-Vij > 0} , and define vectors Wi by Wi = Ui — Y^- Vij. Let Ui denote 
the cardinality of r/j. Now we define M = max{c- Wj : i G /}, S = {i ^ I : c ■ Wi = M} 
and we set a = X]ies'-^i(~l)"*- Note that M simply denotes the highest exponent of 
t appearing in the expansions of the rational functions defined for each i G / in ()1.9p . 
The number a is in fact the sum of the coefficients of t'^^ in these expressions, that is, 
a is the coefficient of t*^ in f{P,t). Now with thes e definitions and notation we can 



state the following result proved by 



Lasserrd (J200J). 



Theorem 1.27. \Lasserre 



200A. Theorem 3.1) 

If c ■ Vij 7^ for alii E I,j E {1, . . . , d}, and if a ^ 0, then M is the optimal value it 
of Integer Program maximize{c ■ x : Ax < b,x > 0,x E 7/}. 



1.5 Summary of results in this thesis 

In this section, we would like to summarize all of the results in this thesis. In Chapter 
121 the first theorem concerns the computation of the toric ideal I a of the matrix A. 

Theorem 1.28. Let A G J/^^'^ and a term order -<w specified by a matrix W . As- 
suming that n and d are fixed, then there are algorithms, that run in polynomial time 
in the size of the input data, to perform the following four tasks: 
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1. Compute a short rational function G which represents the reduced Grobner basis 
of the toric ideal I a with respect to the term order -<]y. 

2. Decide whether the input monomial x" is in normal form with respect to G. 

3. Perform one step of the division algorithm modulo G. 

4- Compute the normal form of the input monomial x"" modulo the Grobner basis 
G. 

The proof of Theorem 11.281 will be given in Section 12.11 S pecial attention will be 



paid to the Projection Theorem (JBarvinok and Wood 



20031 Theorem 1.7) since the 



projection of short rational functions is the most difficult step to implement. Its 
practical efficiency has yet to be investigated. 

Theorem 11.281 can be applied to prove many interesting theorems and corollaries in 
Computational Statistics and Integer Programming. The following corollary can be 
proved by Theorem 11.281 and the fact that a Grobner basis associated to a toric ideal 
I A of the integral matrix A is a Markov basis associate to A. 

Corollary 1.29. Let A G Z'"^'^, where d and n are fixed. There is a polynomial time 
algorithm to compute a multivariate rational generating function for a Markov basis 
M associated to A. This is presented as a short sum of rational functions. 

From Algorithm 11.251 one can see that Theorem II .281 proves the following theorem: 

Corollary 1.30. Let A E Z"^'^, b E Z", and c E Z'^. Given a polyhedron P = 
{x E W^ : Ax = b, X > 0}, compute the mixed integer programming problem, via 
the Grobner basis associated a toric ideal Ia with the term order -<(. obtained by 
Barvinok's rational functions in Theorem \1.2^ 

maximize c ■ x subject to x E P, x E Z*^, 

in polynomial time if we fix n and d. 
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Chapter 121 introduces a new algorithm, improving on Barvinok's original w ork, which 



we ca ll the homogenized Barvinok algorithm. Like the original version in (JBarvinok 



199J), it runs in polynomial time when the dimension is fixed. Then we will apply the 
homogenized Barvinok's algorithm to Commutative Algebra. The Hilbert series of 
S is the rational generating function Xlaes-^"- Barvinok and Woods (2003) showed 
that this Hilbert series can be computed as a short rational generating function in 
polynomial time for fixed dimension. We show that this computation can be done 
without the Projection Theorem (Lemma 12. 4j) when the semigroup is known to be 
normal. 

Theorem 1.31. Under the hypothesis that the ambient dimension d is fixed, 

1. the Ehrhart series of a rational convex polytope given by linear inequalities can 
be computed in polynomial time. The Projection Theorem is not used in the 
algorithm. 

2. The same applies to computing the Hilbert series of a normal semigroup S. 

Chapter |21 discusses Integer and Mixed Integer Programming. We describe a new 
mixed integer programming algorithm via Barvinok's rational functions. This is a 
different approach from the method via the reduced Grobner basis of the toric ideal 
I A associated to the integral matrix A. It is based on Theorem 11.271 in Chapter ^ 
Also in Chapter |21 we give a new algorithm to compute the optimal value and an 
optimal solution for any Integer Linear Program via Barvinok's rational functions. 
In Section IX^ we will show the performance of the BBS algorithm in some knapsack 
problems. Chapter [HI will show a proof of the following theorem: 

Theorem 1.32. Let A G Z™^'', b G Z™, c G Z'^, and assume that number of variables 
d is fixed. Suppose P := {x ^W^ : Ax < 6, x > 0} is a rational convex polytope in 
W^. Given the mixed-integer programming problem 

maximize c ■ x subject to x E {x eM. : x E P, Xi E Tj for i E J C {l,...,(i}}. 
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(A) We can use rational functions to encode the set of vectors (the {A,b,c)-test set): 
{u — V : u is a c — optimal solution, v feasible solution, u,v & Z''}, 

and then solve the MIP problem in time polynomial in the size of the input. 

(B) More strongly, the {A, b, c) test set can be replaced by smaller test sets, such as 
Graver bases or reduced Grobner bases. 

We improve Lasserre's heuristic and give a third deterministic IP algorithm based 
on Barvinok's rational function algorithms, the digging algorithm. In this case the 
algorithm can have an exponential number of steps even for fixed dimension, but 
performs well in practice. See Section IT^ for details. 

Chapter |3] concentrates on computational experiments and explains details of the 
implementation of the software package LattE. We implemented the BBS algorithm 
and the digging algorithm in the second release of the computer software LattE. We 
solved several challenging knapsack problems and compared the performance of LattE 
with the mixed-integer programming solver CPLEX version 6.6. In fact the digging 
algorithm is often surpassed by what we call the single cone digging algorithm. See 
Section HTSl for computational tests. In Section H721 we present some computational ex- 
perience with our current implementation of LattE. We report on experiments with 
families of well-known rational polytopes: multiway contingency tables, knapsack 
type problems, and rational polygons. We demonstrate that LattE competes with 
commercial branch-and-bound software and solves very hard instances and enumer- 
ates some examples that had never been done before. We also tested the performance 
in the case of two-way contingency tables and K ostant's partition function wher e 



special purpose software has been written a 



Beckl(l2nn3[ ): lDe Loera and SturmfeLJ (12001^ 



ready 



Mount 



Baldoni-Silva and Vergnd (|2002l ): 
( 20001 ). In Section lOl we present 



formulas for the Ehrhart quasi-polynomials of several hypersimplices and truncations 
of cubes (e.g. the 24 cell). We show solid evidence that Barvinok's ideas are practi- 
cal and can be used to solve non-trivial problems, both in Integer Programming and 
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Symbolic Computing. 

In Section W^ we present some experimental results with the homogenized Barvinok 
algorithm. It was recently implemented in LattE. Like we will show in Chapter |21 it 
runs in polynomial time when the dimension is fixed. But it performs much better in 
practice (1) when computing the Ehrhart series of polytopes with few facets but many 
vertices; (2) when computing the Hilbert series of normal semigroup rings. We show 
its effectiveness by solving the classical counting problems for 5x5 m,agic squares 
(all row, column and diagonal sums are equal) and 3x3x3x3 magic cubes (all 
line sums in the 4 possible coordinate directions and the sums along main diagonal 
entries are equal). Our computational results are presented in Theorem 14. Ill 
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Chapter 2 



Grobner bases of toric ideals via 
short rational functions 



The main techniques used in this thesis came from the Algebra of polynomial ideals. 
We use special sets of generators called Grobner bases. We deal with ideals associated 
with polyhedra that are called toric ideals (see definitions II. HI . In this chapter we 
present polynomial-time algorithms for computing with toric ideals and semigroup 
rings in fixed dimension. For b ackgrourid on these algebraic objects and their int erplay 



with polyhedral geometry see (jStanlev 



1996 



Sturmfels 



1996; 



Villarreal 



200 iL Our 



results are a direct application of recent results by Barvinok and Woods (2003) on 
short encodings of rational generating functions (such as Hilbert series) . 

2.1 Computing Toric Ideals 



From now on and without loss of generality we will assume that ker{A) fl . 



^>o 



{0}. 



This condition is not restrictive because toric ideal problems can be reduced to this 
particular case via homogenization of the problem. Our assumption implies that for 
all b, the convex polyhedron P = {u & Mf^ : A ■ u = b and m > } is a polytope 
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(i.e. a bounded polytope) or the empty set. We begin by recalling some useful results 
of Barvinok and Woods (2003): 



Lemma 2.1. \Barvinok and WoodL \2003i. Theorem 3.6) Let 81,82 be finite subsets 
of 'Z'^, for d fixed. Let f{8i,x) and f{82,x) be their generating functions, given as 
short rational functions with at most k binomials in each denominator. Then there 
exists a polynomial time algorithm, which, given f{8i,x), computes 



fi8,n82,x) = E^^'n 



(1 -x^")---(l -x^'O 
with s < 2k, where the 7j are rational numbers, Ui,Vij are nonzero integer vectors, 
and I is a polynomial- size index set. 

The following lemma was proved by Barvinok and Woods using Lemma f2. 11 



Lemma 2.2. Warvinok and Woods . 120031. Corollary 3. 7) Let 81, 82, ■ ■ ■ ,8m be finite 
subsets of Z'^, for d fixed. Let f{8i, x) for i = 1 . . .m be their generating functions, 
given as short rational functions with at most k binomials in each denominator. Then 
there exists a polynomial time algorithm, in the input size, which computes 



f{8iU82U...8m,x) = $^7..-7Y 



with s < 2k, where the 7^ are rational numbers, Ui, Vij are nonzero integer vectors, 
and I is a polynomial- size index set. Similarly one can compute in polynomial time 
f{8i\82,x) as a short rational function. 

We will use the Intersection Lemma and the Boolean Operation Lemma to extract 
special monomials present in the expansion of a generating function. The essential 
step in the intersection algorithm is the use of the Hadamard product (see Algorithm 
II. 3|) and a special monomial substitution. The Hadamard product is a bilinear op- 
eration on rational functions (we denote it by *). The computation is carried out 
for pairs of summands as in ()1.2|) . Note that the Hadamard product mi * m2 of two 
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monomials ?Tii,m2 is zero unless mi = 7712- We present an example of computing 
intersections. 

Example 2.3. Let S'j = { x G M : z - 2 < x < i } fl Z for z = 1,2. We rewrite their 



rational generating functions as in the proof of Theorem 3.6 in (JBarvinok and Woods . 



m^: f{Si,z) = Tf^ + TT^f-TT = 71^ + 7T3f-TT = 9n + ^12, and f{S2,z) 



{l-z) ^ (1-z-i) ~ (1-2-1) -r (i_2-i) 



1 z ~z^ z 

(T^ + (1-2-1) = (i„^-i) + (i_^-i) = fi'21 + fi'22- 

We need to compute four Hadamard products between rational functions (^jj, whose de- 
nomi nators are products of binomials and whose numerators are monomials. Lemma 



3.4 in lBarvinok and Woodsl ()2nn3( ) says that, these Hadamard products are essentially 
the same as computing the rational function, as in Equation ()1.2|1 . of the auxiliary 
polyhedron {(ei,e2)|pi -|- a^e^ = P2 + «2e27 Q > 0}. Here Pi,P2 are the exponents of 
numerators of gij 's involved and ai, 02 are the exponents of the binomial denomina- 
tors. For example, the Hadamard product gu * (722 corresponds to the polyhedron 

— 2 

{(ei, e2)|e2 = 4 + ei, e^ > 0}. The contribution of this half line is —7^37^- We find 

f{Sl,z)* f{S2,z) = 911*921+912*922 + 912*921+911*922 



2-2 


2-1 2-2 


(1-2-1) 1 (1_^-1) 


(1-2-1) (1-2-1) 


z-z-'^ _ 1 1 ^ 

1-2-1 - l + Z 


= f{SinS2,z) 



Another key subroutine introduced by Barvinok and Woods is the following Projection 
Theorem. In Lemmas 12.11 12.21 and 12. 4t the dimension d is assumed to be fixed. 



Lemma 2.4. Warvinok and Woodi . 



2003 . Theorem 1.7) Assume the dimension d is 



a fixed constant. Consider a rational polytope P C M'^ and a linear map T : Z"' ^ Z'^ . 
There is a polynomial time algorithm which computes a short representation of the 
generating function / (T(P fl Z'^), x) . 

Defining a term order -<w hj a. d x d integral matrix W (see details in Section [OJ, 
we have the following lemma. 



CHAPTER 2. 



Grobner bases of toric ideals via short rational functions 



30 



Lemma 2.5. Let S C Z^ be a finite set of lattice points in the positive orthant. 
Suppose the polynomial f{S,x) = '^s^gX^ is represented as a short rational function 
and let -<vi/ be a term order. We can extract the (unique) leading monomial of f{S,x) 
with respect to -<w in polynomial time. 

Proof: The term order -<y/ is represented by an integer matrix W . For each of the rows 
Wj of W we perform a monomial substitution xi := x'jV"^\ Note that t is a "dummy 
variable" that we will use to keep track of elimination. Such a inononi ial substitution 
can be computed in polynomial time by (JBarvinok and WoodsL 120031 Theorem 2.6). 
The effect is that the polynomial f{S,x) gets replaced by a polynomial in the t and 
the x's. After each substitution we determine the degree in t. This is done as follows: 
We want to do calculations in univariate polynomials since this is faster so we consider 
the polynomial g{t) = f{S, l,t), where all variables except t are set to the constant 
one. Clearly the degree of g{t) in t is the same as the degree of f{S,x',t). We 
create the interval polynomial i[p^q]{t) = Yl'i=p^^ which obviously has a short rational 
function representation. Compute the Hadamard product of i[p.q](t) with g(t). This 
yields those monomials whose degree in the variable t lies between p and q. We will 
keep shrinking the interval [p, q] until we find the degree. We need a bound for the 
degree in t of g{t) to start a binary search. An upper bo und U can be found via 



linear programming or via the estimate in Theorem 3.1 of ()Lasserrd . |200J) which is 
an easy manipulation of the numerator and denominator of the fractions in g{t). It 
is clear that log(f/) is polynomially bounded. In no more than log{U) steps one can 
determine the degree in t of f{S,x,t) by using a standard binary search algorithm. 
Let a be a polynomial-size upper bound on the highest total degree of a monomial 
appearing in the genera ting func t ion f (S,x). We can again apply linear program- 



ming or the estimate of (JLasserre 



20041 ) to compute such an a (just as we computed 



U before). Once the highest degree r in t is known, we compute the Hadamard 
product of f{S,x,t) and t^/i(x), where h{x) is the rational generating function en- 
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coding the lattice points contained inside the box [0, of. This will capture only 
the desired monomials. Then compute the limit as t approaches 1. This can be 
done in polynomial time using residue techniques. The limit represents the subseries 
H{S,x) = J2i3-w=r^^- Repeat the monomial and highest degree search for the row 
Wj+i,Wj+2, etc. Since -<w is a term order, after doing this d times we will have only 
one single monomial left, the desired leading monomial. D 

One has to be careful when using earlier Lemmas (especially the projection theorem) 
that the sets in question are finite. We need the following well-known bound: 



Lemma 2.6. nSturmfelA . \l99(\ . Lemma 4-6 and Theorem, J^.l) Let M he equal to 
{n + l){d — n)D{A), where A is an n x d integral matrix and D{A) is the biggest 
n X n suhdeterminant of A in absolute value. Any entry of an exponent vector of any 
reduced Grobner basis for the toric ideal I a is less than M . 

Proposition 2.7. Let A G Z"^"^, W G TJ^^'^ specifying a term order -<w Assume 
that n and d are fixed. 

1) There is a polynomial time algorithm to compute a short rational function G which 
represents a universal Grobner basis of Ia- 

2) Suppose we are given the term order -<y/ and a short rational function encoding 
a finite set of binomials x"^ — x"" now expressed as the sum of monomials Y^x^y"". 
Assume M is an integer positive bound on the degree of any variable for any of the 
monomials. One can compute in polynomial time a short rational function encoding 
only those binomials x" — x"" that satisfy x^ -<w x^ ■ 

3) Suppose we are given a sum of short rational functions f{x) which is identical, in 
its monomial expansion, to a single monomial x"" . Then in polynomial time we can 
recover the (unique) exponent vector a. 

Proof: 1) Set M = {n + l){d — n)D{A) where D{A) is again the larges t absolute 



value of any n x n-sub determinant of A. Using Barvinok's algorithm in (|Barvinok 
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1994? ) ■ we compute the following generating function in 2d variables: 

G{x,y) = ^{xV : Au = Av andO<Ui,Vi<M}. 

This is the sum over all lattice points in a rational polytope. Lemma 171)1 above implies 
that the toric ideal I a is generated by the finite set of binomials x^ — x^ corresponding 
to the terms x'^y'" in G{x, y). Moreover, these binomials are a universal Grobner basis 

oUa. 

2) Denote by Wi the i-th row of the matrix W which specifies the term order. Suppose 
we are given a short rational generating function Go{x,y) = ^x'^y^' representing 
a set of binomials x" — x^ in J^, for instance Gq = G in part (1). In the following 
steps, we will alter the series so that a term x'^y'" gets removed whenever u is not 
bigger than v in the term order -<w/. Starting with Hq = Gq, we perform Hadamard 
products with short rational functions f{S] x, y) for S C Z^'^. 

Set Hi = Hi^i * /({(m, f) : WiU = WiV, < Uj,Vj < M, j = l...d}), and Gi = 
Hi-i * f{{{u,v) : WiU > WiV + 1, < Uj^Vj < M j = l...d}). All monomials 
x'^y'" G Gj have the property that WiU = WiV for i < j, WjU > WjV, and thus v ^w u. 
On the other hand, if v -<w u then there is some j such that WiU = WiV for i < j, 
WjU > WjV, and we can conclude that x^y^ G Gj. Note that if = Gi U G2 U . . . U G^ 
is actually a disjoint union of sets. The rational function that gives the union, can 
be computed in polynomial time by Lemma [2.21 In practice, the rational generating 
functions representing the Gj's can be simply added together. The short rational 
function H encodes exactly those binomials in Go that are correctly ordered with 
respect to ^ly- We have proved our claim since all of the above constructions can be 
done in polynomial time. 

3) Given f{x) we can compute in polynomial time the partial derivative df{x)/dxi. 
This puts the exponent of Xj as a coefficient of the unique monomial. Computing 
the derivative can be done in polynomial time by the quotient and product derivative 
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rules. Each time we differentiate a short rational function of the form 



j,bi 



(1 - X=i'')(l - X"^-') ... (1 - X'^'i-^) 

we add polynomially many (binomial type) factors to the numerator. The factors 
in the numerators should be expanded into monomials to have again summands in 
short rational canonical form 7- — ctttt, — %' ^ n ,-, — ctt-t- Note that at most 2*^ monomials 
appear each time {d is a constant). Finally, if we take the limit when all variables Xi 
go to one we will get the desired exponent. D 

Example 2.8. Using LattE we compute the set of all binomials of degree less than or 

r 1 1 1 1 1 
equal 10000 in the toric ideal I a of the matrix A = . This matrix rep- 

[ 1 2 3 J 

resents the Twisted Cubic Curve in algebraic geometry. We find that there are exactly 

195281738790588958143425 such binomials. Each binomial is encoded as a monomial 
x\^X2^xl^x^*y\^y2^y^^y\*. The computation takes about 40 seconds. The output is a 
sum of 538 simple rational functions of the form a monomial divided by a product 
such as (1 - sm) (1 - 2^) (1 - xm) (1 - x^xm^) (1 - 2:3^3) (1 - X2y2). □ 

Proof of Theorem ITT^ 

The proof of Theorem ll. 281 will require us to project and intersect sets of lattice points 
represented by rational functions. We cannot, in principle, do those operations for 
infinite sets of lattice points. Fortunately, in our setting it is possible to restrict our 
attention to finite sets. Besides Lemma IT6l for the size of exponents of Grobner bases, 
we need a bound for the exponents of normal form monomials: 

Lemma 2.9. Let x" he the normal form of x"" with respect to the reduced Crobner 
basis G of a toric ideal I a for the term order -<y/ (associated to the matrix W). Every 
coordinate of u is bounded above by L = {n + l)dD{A)d, where D{A) is the biggest 
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subdeterminant of A in absolute value, a denotes the largest coordinate of the exponent 
vector a. 

Proof: We note that u is a point in the (bounded) convex polytope defined by the 
following inequalities in v. Av = Aa, and f > (it is forced to be bounded for all 
a because we assumed ker{A) n M>o = {0}). Thus each coordinate of u is bounded 
above by the corresponding coordinate of some vertex of this polytope. Let v be 
such a vertex. The non-zero entries of v are given by B~^Aa where i? is a maximal 
non-singular square submatrix of A. Clearly, each entry of B^^A is bounded above 
by D{A), and hence each entry of v is bounded above by L. We conclude that L is 
an upper bound for the coordinates of u. U 

Proof of Theorem \1.2t^ Proposition 12.71 gives a Grobner basis for the toric ideal I a 
in polynomial time. We now show how to get the reduced Grobner basis from it in 
three easy polynomial time steps. The input is the the n x d integral matrix A and 
the d X d term order matrix W. The algorithm for claim (1) of Theorem 11.281 has 
three steps: 

Step 1. Let M be equal to {n + l){d — n)D{A), as in Lemma 12.61 for given input 
matrix A. As in Proposition 12.71 compute the generating function which encodes 
binomials of highest degree M on variables that generate /^: 

/(x, y) = ^{ x^'y'' : Au = Av and < Uj, Vj < M ioi j = 1 . . .d}, 

Next we wish to remove from f{x, y) all incorrectly ordered binomials (i.e. those 
monomials x'^y^ with u -<w v instead of the other way around). We do this using 
part 2 of Proposition 12.71 We obtain from it a collection Gq,Gi, . . . ,Gd of rational 
functions encoding disjoint sets of lattice points. We call /(x, y) the generating func- 
tion representing the union of Gq, . . . , Gd- This can be computed in polynomial time 
by adding the rational functions of the Gi together (since they are disjoint). The 
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reader should notice that this updated f{x, y) contains only those monomials of the 

old f{x,y) that are now correctly ordered. 

Let gi{x) be the projection of Gi onto the first group of x- variables and denote by g{x) 

the rational function that represents the union of the gi{x). The rational function g{x) 

can be computed in polynomial time by the projection theorem of Barvinok-Woods, 

i.e. Lemma 12.41 It is important to note that g{x) is the result of projecting f{x,y) 

into the first group of variables. This is true because a linear projection of the union 

of disjoint lattice point sets (i.e. those represented by Gi) equals the union of the 

projections of the individual sets. In conclusion, g{x) is the sum over all non-standard 

monomials having degree at most M in any variable. 

Step 2. Write r{x,M) = n(i3^ + '^' -i ) for the generating function of all x- 

i=l * 

monomials having degree at most M in any variable. Note that this is a large, but 
finite, set of monomials. We compute the following Hadamard product of d rational 
functions in x and Boolean complements (we denote them by \): 

r(x, M)\xi ■ g{x) 1 * I r(x, M)\x2 ■ g{x) J * ■ ■ ■ * I r(x, M)\xd ■ g{x) J . 

This is the generating function over those monomials all of whose proper factors are 
standard modulo the toric ideal I a and whose degree in any variable is at most M. 
Step 3. Let h{x,y) denote the ordinary product of the resulting rational function 
from Step 2 with 

r{y, M)\g{y) = /_^{ y^ '■ v standard monomial modulo /a of highest degree M^. 

Thus h{x, y) is the sum of all monomials x'^y^ such that x^ is standard and x" is a 
monomial all of whose proper factors are standard monomials modulo the toric ideal 
I A and, finally, the highest degree in any variable is at most M. 
Compute the Hadamard product G{x, y) := /(x, y) * h{x, y). This is a short rational 
representation of a polynomial, namely, it is the sum over all monomials x'^y'" such 
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that the binomial x" — x'" is in the reduced Grobner basis of Ia with respect to W 
and x^ -<w x"^- This completes the proof of the first claim of Theorem 11.281 
We next give the algorithm that solves claims 2 and 3 of Theorem 11.281 This will be 
done in four steps (1,2,3,4). We are given an input monomial x" for which we aim to 
determine whether it is already in normal form. 

Step 4 Perform Steps 1,2,3. Let G{x,y) be the reduced Grobner basis of Ia with 
respect to the term order W encoded by the rational function obtained at the end of 
Step 3. Let r(x, a) be, as before, the rational function of all monomials having degree 
less than a on any variable. Thus G'{x, y) = r{x, a) ■ G{x, y) consists of all monomials 
of the form x^{x^y'^') where x" — x^' is a binomial of the Grobner basis and where 
< s < a. Thus x'^x^ is a monomial divisible by some leading term of the Grobner 
basis. 

Given a monomial x"' consider b{x,y), the rational function representing the lattice 
points of {{u,v) : u = a, < Vj < L for j = l...d}. The Hadamard product 
G{x, y) = G'{x, y) * b{x, y) is computable in polynomial time and corresponds to 
those binomials in G{x,y) that can reduce x"'. If G{x,y) is empty then x" is in 
normal form already, otherwise we use Lemma 12.51 and part 3 of Proposition 12.71 to 
find an element x^y^ G G{x,y) and reduce x" to x'^~"+'". We may assume that the 
coefficient of the encoded monomial is one, because we can compute the coefficient in 
polynomial time using residue techniques, and divide our rational function through 
by it. 

Finally, we present the algorithm for claim 4 in Theorem ll.28l in four steps (1,2,5,6). A 
curious byproduct of representing Grobner bases with short rational functions is that 
the reduction to normal form need not be done by dividing several times anymore. 
Step 5. Redo all the calculations of the Steps 1,2,3 using L = [n + l)dD{A)a 
from Lemma 12.91 instead of M. Note that the logarithm of L is still bounded by a 
polynomial in the size of the input data {A, W, a). Let f{x, y) and g{x) from Step 1,2 



CHAPTER 2. Grobner bases of toric ideals via short rational functions 37 

(now recomputed with the new bound L) and compute the Hadamard product 

H{x,y) := f{x,y)* (r{x,L) ■ [r{y,L)\g{y)) 

This is the sum over all monomials x'^y'" where x"" is the normal form of x" and highest 
degree of x" on any variable is L. Since we took a high enough degree, by Lemma 
12.91 the monomial x"'y''^, with x'^ the normal form of x"', is sure to be present. 
Step 6. We use H{x, y) as one would use a traditional Grobner basis of the ideal I^. 
The normal form of a monomial x"' is computed by forming the Hadamard product 
H{x,y) * (x" ■ r{y,L)). Since this is strictly speaking a sum of rational functions 
equal to a single monomial, applying Part 3 of Proposition 12.71 completes the proof 
of Theorem ir^ D 

2.2 Computing Normal Semigroup Rings 

We will show in Chapter ID th at a major practical bottleneck of the original Barvinok 



algorithm in (JBarvinok 



1994 ) is the fact that a polytope may have too many vertices. 
Since originally one visits each vertex to compute a rational function at each tangent 
cone, the result can be costly. For example, the well-known polytope of semi-magic 
cubes in the 4x4x4 case has over two million vertices, but only 64 linear inequalities 
describe the polytope. In such cases we propose a homogenization of Barvinok's 
algorithm working with a single cone. 
here is a second motivation for looking at the homogenization. Barvinok and Woods 



tearvinok and Woods . 



20031 ) proved that the Hilbert series of semigroup rings can be 
computed in polynomial time. We show that for normal semigroup rings this can be 
done simpler and more directly, without using the Projection Theorem. 
Given a rational polytope P in M'^"^, we set i{P,m) = ^{z G Z''"^ : z G mP}. The 
Ehrhart series of P is the generating function ^^=q ^(-P, m)t"^. 
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Algorithm 2.10 (Homogenized Barvinok algorithm). 

Input: A full- dimensional, rational convex polytope P in M.'^~^ specified by linear 

inequalities and linear equations. 

Output: The Ehrhart series of P. 

1. Place the polytope P into the hyperplane defined by x^ = 1 in M'^. Let K be 
the (i-dimensional cone over P, that is, K = cone{{{p, 1) : j» G P})- 

2. Compute the polar cone K*. The normal vectors of the facets of K are exactly 
the extreme rays of K*. If the polytope P has far fewer facets then vertices, 
then the number of rays of the cone K* is small. 

3. Apply Barvinok's decomposition of K* into unimo dular cones. Polarize back 
each of these cones. It is known, e.g. Corollary 2.8 in f|Barvinok and Pommersheim 



19991 ). that by dualizing back we get a unimodular cone decomposition of K. All 
these cones have the same dimension as K. Retrieve a signed sum of multivariate 
rational functions which represents the series J^aeKnZ'' ^"• 

4. Replace the variables Xj by 1 for z < rf — 1 and output the resulting series in 
t = Xd- This can be done using the methods in Chapter |3 

One of the key steps in Barvinok's algorithm is that any cone can be decomposed as 
the signed sum of (indicator functions of) unimodular cones. We will talk about this 
in detail on Section 14.1.11 Chapter HI 



Theorem 2.11 (see (JBarvinokl . Il994l )^. Fix the dimension d. Then there exists a 
polynomial time algorithm which decomposes a rational polyhedral cone /T C M^ into 
unimodular cones Ki with numbers ej G { — 1,1} such that 

f{K) = Y,^.J{Ki), |/|<cx). 

i&I 

Originally, Barvinok had pasted together such formulas, one for each vertex of a 
polytope, using a result of Brion. Using Algorithm 12. 101 we can prove Theorem 11.311 
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Proof of Theorem ] 1.!^ 11 We first prove part (1). The algorithm solving the problems 
is Algorithm 12. lUI Steps 1 and 2 are polynomial when the dimension is fixed. Step 
3 follows from Theorem 12.111 We require a special monomial substitution, possibly 



with some poles. This can be done in polynomial time by (JBarvinok and Woods . 



2mm. 

Part (2): Recall the characterization of the integral closure of the semigroup S as the 
intersection of a pointed polyhedral cone with the lattice Z'^. From this it is clear that 
Algorithm I2.1UI computes the desired Hilbert series, with the only modification that 
the input cone K is given by the rays of the cone and not the facet inequalities. The 
rays are the generators of the monomial algebra. But, in fixed dimension, one can 
transfer from the extreme rays representation of the cone to the facet representation 
of the cone in polynomial time. D 
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Chapter 3 

Theoretical applications of rational 
functions to Mixed Integer 
Programming 

We now discuss how all these ideas can be used in Discrete Optimization. 

3.1 The {A,b,c) test set algorithm 

In all our discussions below, the input data are a.n m x d integral matrix A and 
an integral ?n-vector b. For simplicity we assume it describes a polytope P = {x ^ 
M.'^\Ax < b,x > 0}. We assume that there are no redundant inequalities and no hidden 
equations in the system. This polytope P = {x ^ M.'^\Ax < fe, a; > 0} is equivalent 
to the expression of {x G M'^jAx = b,x > 0}, where A G Z™^*^ and d = m + d. 
We can transform the expression {x G M'^lAx = 6, x > 0} to the expression {x G 
M'^lAx < b,x > 0} hj projecting down P to a full dimensional polytope with Hermite 
normal form and we can also transform the expression {x G M'^jAx < 6, x > 0} to the 
expression {x G M"'|ylx = b,x > 0} hj introducing slack variables. 



CHAPTER 3. 



Theoretical applications of rational functions to MIP 



41 



First we would like to remind a reader of Barvinok's rational functions (see details on 
Chapter nj. By Theorem ll.il with a given P, if we fix d there is a polynomial time 
algorithm to compute Barvinok's short rational functions in the form of 



f{P,z) = Y,E^ 



(3.i: 






where / is a polynomial sized finite indexing set, and where Ei G {1,-1} and Uj, Vij G 
1/ for all i and j. In this section we will show how to apply Barvinok's short rational 
functions in ()3.1|) to Mixed Integer Programming. 

Proof of Theorem \1.3^ - We only show the proof of part (A). The proof of part (B) 
appears in Chapter |21 We first explain how to solve integer programs (where all 
variables are demanded to be integral). This p a rt of t he proof is essentially the proof 



of Lemma 3.1 given in 



Hosten and Sturmfelsl (J20031 ) for the case Ax = b, x > 0, 



instead of Ax < b, but we emphasize the fact that b is fixed here. We will see how the 
techniques can be extended to mixed integer programs later. For a positive integer 
R, let 

r{x,R)=l[i^ 



tXj A 



i=l 



X, 



Xi 1 

be the generating function encoding all x-monomials in the positive orthant, having 
degree at most R in any variable. Note that this is a large, but finite, set of monomials. 
Suppose P = ^x E M,'^ : Ax < b, x > 0} i s a non empty polytope. Using Barvinok's 



algorithm in 



Barvinok and PommersheimI (jl999f ). compute the following generating 



function in 2d variables: 

f{x, y) = ^{ x^y" : Au < b, Av < b, u, v > 0, c ■ u - c ■ v > 1, and u,v e Z'^}. 

This is possible because we are clearly dealing with the lattice points of a rational 
polytope. The monomial expansion of f{x,y) exhibits a clear order on the variables: 
x^y" where c-u > c-v. Hence v is not an optimal solution. In fact, optimal solutions 
will never appear as exponents in the y variables. 
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Now let g{y) be the projection of f{x,y) onto the y-variables variables. Thus g{y) 
is encoding all non-optimal feasible integral vectors (because the exponent vectors of 
the x's are better feasible solutions, by construction), and it can be computed from 
f{x, y) in polynomial time by Lemma EU Let V{P) be the vertex set of P and choose 
an integer R > max{fj : v G V{P), 1 < i < d} (we can find such an integer R via 
linear programming). Define f{x,y) and g{x) as above and compute the Hadamard 
product 

Hix,y) := f{x,y)*[{r{x,R)- g{x))r{y,R)]. 

This is the sum over all monomials x'^y'" where u, v E P and where u is an optimal 
solution. The reader should note that the vectors u — v form a test set (an enormous 
one), since they can be used to improve from any feasible non-optimal solution v. 
This set is what we called the (A, 6, c)-test set. It should be noted that one may 
replace H{x, y) by a similar encoding of other test sets, like the Graver test set or a 
Grobner basis (see Chapter |^ for details). 

We now use H{x, y) as one would use a traditional test set for finding an optimal 
solution: Find a feasible solution a inside the polytope P using Lemma 12.51 and 
Barvinok's Equation ()3.1|) . Improve or augment to an optimal solution by computing 
the Hadamard product 

H{x,y)*{y\{x,R)). 

The result is the set of monomials of the form x"?/" where u is an optimal solution. 
One monomial of the set, say the lexicographic largest, can be obtained by applying 
Lemma 12.51 This concludes the proof of the case when all variables are integral. 
Now we look at the mixed integer programming case, where only Xi with i G J C 
{1, . . . , d} are required to be integral. Without loss of generality, we may assume that 
J = {r, . . . , (i} for some r,l <r < d. Thus, splitting A into {B\C), we may write the 
polytope P as {(x,x') : Bx + Cx' < b, x,x' > 0} where the variables corresponding 
to B are not demanded to be integral. Consider a vertex optimal solution x to the 
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mixed integer problem. The first key observation is that its fractional part can be 
written as xj = B^^ih — Cx') where h — Cx' is an integer vector. Here B~^ denotes 
the inverse of a submatrix of B. This follows from the theory of linear programming, 
when we solve the mixed integer program for fixed x' = x! . 

The denominators appearing are then contributed by B~^ . Then every appearing 
denominator is a factor of M, the least common multiple of all determinants of a 
square submatrix of A. It is clear M can be computed in polynomial time in the 
size of the input. This complexity bound holds, since the number of such square 
submatrices is bounded by a polynomial in m, the number of rows of A, of degree d, 
the number of columns of A. Moreover, each of these determinants can be computed 
in time polynomial in the size of the input, and therefore, M itself can be computed 
in time polynomial in the size of the input in fixed dimension d. Thanks to this 
information, we know that if we dilate the original polytope P by M , the optimal 
solutions of the mixed integer program become, in the dilation MP ^ optimal integral 
solutions of the problem 

maximize c ■ x subject to x G MP^ x G Z 

with the additional condition that the coordinates with index in J are multiples of 
M . Ignoring this condition involving multiples of M for a moment, we see that, as we 
did before, we can obtain an encoding of all optimal improvements as a generating 
function H{x,y). 



Let r{x, R) 






To extract those vec- 



tors whose coordinates indexed by J are multiples of M, we only need to intersect 
(Hadamard product again) our generating function H{x, y) with the generating func- 
tion f{x,R)f{y,R). Then only those vectors whose coordinates indexed by J are 
multiples of M remain. This completes the proof of the theorem. D 
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3.2 The Digging Algorithm 

In what follows we present a strengthening of Lasserre's heuristic and discuss how 
to use Barvinok's short rational functions to solve integer programs using digging. 
Suppose A G Z™^"^, b G Z™ and finite S C Z'^ are given. We consider the family of 
integer programming problems of the form maximize {c ■ x : Ax < b,x > 0,x G Z*^}, 
where c ^S. We assume that the input system of inequalities Ax < b,x > defines 
a polytope P C W^, such that P fl Z'^ is nonempty. 

When the hypotheses of Theorem 11.271 are met, from an easy inspection, we could 
recover the optimal value of an integer program. If we assume that c is chosen 
randomly from some large cube in Z'', then the first condition is easy to obtain. 
Unfortunately, our computational experiments (see Section 14. 5|) indicate that the 
condition cr 7^ is satisfied only occasionally. Thus an improvement on the approach 
that Lasserre proposed is needed to make the heuristic terminate in all instances. Here 
we explain the details of an algorithm that digs for the coefficient of the next highest 
appearing exponent of t. For simplicity our explanation assumes the easy-to-achieve 
condition c ■ Vij ^ 0, for all Vij. 

As before, take Equation (j3.ip computed via Barvinok's algorithm. Now, for the 
given c, we make the substitutions Zk = Vkt'^'', for k = 1,. . . ,d. Then substitution 
into (|3.ip yields a sum of multivariate rational functions in the vector variable y and 
scalar variable t: 

g{P; y,t) = J2 ^^FFd— T :■ (3-2) 

On the other hand, the substitution on the left-side of Equation ()3.1|1 gives the fol- 
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lowing sum of monomials. 

g{P-,y.t)= J2 ^"^""- (3-3) 

aePnzd 

Both equations, (jH.Hj) and ()3.2|1 . represent the same function g{P;y,t). Thus, if we 
compute a Laurent series of (j3.2j) that shares a region of convergence with the series in 
(13 .3^ . then the corresponding coefficients of both series must be equal. In particular, 
because P is a polytope, the series in (j3.3|) converges almost everywhere. Thus if we 
compute a Laurent series of (j3.2j) that has any nonempty region of convergence, then 
the corresponding coefficients of both series must be equal. Barvinok's algorithm 
provides us with the right hand side of (j3.2|) . We need to obtain the coefficient of 
highest degree in t from the expanded Equation (refeq:d). We compute a Laurent 
series for it using the following procedure: Apply the identity 



1 —y-'"ijf-C-Vij 

^ (3.4) 



to Equation ()3.2|1 . so that any Vij such that c ■ Vij > can be changed in "sign" to be 
sure that, for all Vij in ()3.2j) . c ■ Vij < is satisfied (we may have to change some of 
the Ei, Ui and Vy using our identity, but we abuse notation and still refer to the new 
signs as Ei and the new numerator vectors as Ui and the new denominator vectors as 
Vij). Then, for each of the rational functions in the sum of Equation ()3.2|) compute a 
Laurent series of the form 

d 
Ei y^'f"' JJ(1 + y^^H""-"^^ + {y^'if^'if + ...). (3.5) 

i=i 
Multiply out each such product of series and add the resultant series. This yields 

precisely the Laurent series in ()3.3|) . Thus, we have an algorithm to solve integer 

programs: 

Algorithm: {Digging Algorithm): 
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Input: A e Z"^^'^, 6 G Z"^, c G S. 

Output: optimal value and optimal solution of maximize{c ■ x : Ax < b,x > 0,x & 

U^} for all c^^. 

Procedure: for each c G H, do 

1. Use the identity ()3.4|) as necessary to enforce that all Vij in ()3.2|) satisfy c-Vij < 0. 

2. Via the expansion formulas ()3.5p . find ()3.3|) by calculating the terms' coefficients. 
Proceed in decreasing order with respect to the degree of t. This can be done 
because, for each series appearing in the expansion formulas fl3.5|) . the terms of 
the series are given in decreasing order with respect to the degree of t. 

3. Continue calculating the terms of the expansion ()3.3|) . in decreasing order with 
respect to the degree of t, until a degree n of t is found such that for some a G Z'', 
the coefficient of y^t" is non-zero in the expansion ()3.3p . 

4. Return "n" as the optimal value of the integer program and return a as an 
optimal solution. 

We close this section by noticing that one nice feature of the digging algorithm is 
if one needs to solve a family of integer programs where only the cost vector c is 
changing, then Equation (j3.2|) can be computed once and then apply the steps of the 
algorithm above for each cost vector to obtain all the optimal values. 
Given the polytope P := {x & M.'^ : Ax < b, x > 0}, the tangent cone K^ at a 
vertex v of P is the pointed polyhedral cone defined by the inequalities of P that 
turn into equalities when evaluated at v. We will show in Chapter ID th at a major 



Barvinokl (|199J) is the fact 



practical bottleneck of the original Barvinok algorithm in 
that a polytope may have too many vertices. Since originally one visits each vertex to 
compute a rational function at each tangent cone, the result can be costly. Therefore 
a natural idea for improving the digging algorithm is to compute with a single tangent 



CHAPTER 3. Theoretical applications of rational functions to MIP 47 

cone of the polytope and revisit the above calculation for a smaller sum of rational 
functions. Let vertex v* give the optimal value for the given linear programming 
problem and we only deal with the tangent cone K^* . Suppose we have the following 
integer programming problem: 

(IP) maximize c ■ x subject to x G P fl Z*^, 

where P := {x eW^ : Ax <h, x >Q], A e Z"*^^ and h G Z"*. 

Then we have the following linear programming relaxation problem for the given 

integer programming problem: 

(LP) maximize c ■ x subject to x E P. 



One of the vertices of P gives the optimal value for (LP) ISchriiverl (|1986I ). Let V{P) 

be the vertex set of P and v G V{P) be a vertex such that c ■ i; is the optimal value 

for (LP). Then, clearly, the tangent cone Ky at v contains P. So, if we can find an 

integral point x* G K^ such that c ■ x* > c ■ x, "ix E P and x* E P Cl Z"', then x* is 

an optimal solution for (IP). The outline for the single cone digging algorithm is the 

following: 

Algorithm: {Single Cone Digging Algorithm): 

Input: A G Z™^"', b G Z™, c G Z^. 

Output: optimal value and optimal solution of maximize{c ■ x : Ax < 6, x > 0,x G 

In the following steps, we replace P by Ky in the notation. 

1. Compute a vertex f of P such that c ■ v = maximize{c ■ x : Ax < b,x > 0}. 

2. Compute the tangent cone Ky at v and compute the short rational function (j3.2j) 
encoding the lattice points inside K^. 

3. Use the identity ()3.4|) as necessary to enforce that all Vij in ()3.2|) satisfy c-Vij < 0. 
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4. Via the expansion formulas (j3.5|) . find (J3.3J1 by calculating the terms' coefficients. 
Proceed in decreasing order with respect to the degree of t. This can be done 
because, for each series appearing in the expansion formulas (13 .51) . the terms of 
the series are given in decreasing order with respect to the degree of t. 

5. Continue calculating the terms of the expansion ()3.3|) . in decreasing order with 
respect to the degree of t, until a degree ra of t is found such that: 

• for some a G Z'^, the coefficient of ?/"t" is non-zero in the expansion ()3.3|) . 

• Aa <b, a>0. 

6. Return "n" as the optimal value of the integer program and return a as an 
optimal solution. 

From Table 14.171 and Table 14.181 one can find that the single cone digging algorithm 
is very practical compared to the BBS algorithm and the original digging algorithm. 
This algorithm is faster and more memory efficient than the original digging algorithm 
in practice, since the number of unimodular cones for the single cone digging algorithm 
is much less than the number of unimodular cones for the original digging algorithm. 
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Chapter 4 



Experimental results: development 

of LattE 



4.1 LattE's implementation of Barvinok's algorithm 



In this section, we go through the steps of Barvinok's algorithm, showing how we 
implemented them in LattE. Barvinok's algorithm relies on two important new ideas: 
the use of rational functions as efficient data structures and the signed decompositions 
of cones into unimodular cones. 

Let P C M'^ be a rational convex polyhedron and let f{P, z) be the multivariate 
generating function defined in (jl.lj) . Let f be a vertex of P. Then, the supporting cone 
K{P, v) of P at f is K{P, v) = v + {u eW'- : v + 5u & P ioi all sufficiently small 5 > 
0}. Let V{P) be the vertex set of P. One crucial component of Barvinok's algorithm 
is the ability to distribute the computatio n on the ver tices of the polytope. This 



follows from the seminal theorem of Brion ( Brion 



im): 



Theorem 4.1. Wrion 



1983} Let P be a rational polyhedra and let V{P) he the vertex 
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set of P. Then, 



f{P,z)= J2 fiKiP,v),z). 

v£V(P) 



Example 4.2. Consider the integral quadrilateral shown in Figure \4^ The vertices 
are Vi = (0,0), V2 = (5,0), V3 = (4,2), andV^ = (0,2). 
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Figure 4.1: A quadrilateral in Example 14. 21 



We obtain four rational generation functions whose formulas are 

f{Kv,,z) = ,^[':^XT!::%v /(^v.,^) 



(l-22-l){l-2l)' 



Indeed, the result of adding the rational functions is equal to the polynomial 

Zi^ + Zi'^Z2 + Zi'^ + 2:1^2:2^ + 222:1^ + Zi^ + Zi^Z2^ + Z2Zi'^ + Zi^ + Zi^ Z2^ + Z1Z2 



Z\ 



ZiZ2^ + Z2^ + 22 



1. 



D 



In order to use Brion's theorem for counting lattice points in convex polyhedra, we 
need to know how to compute the rational generating function of convex rational 
pointed cones. For polyhedral cones this generating function is a rational function 



whose numerator and denominator have a we 
in 



StanlevI ()l997l Chapter 4) and in 



^understood geometric meaning (see 



StanlevI ()1980l Corollary 4.6.8) for a clear ex- 



planation). We already have a "simple" formula when the cone is a simple cone: 
Let {mi,-U2, • • • 5 Wfc} be a set of linearly independent integral vectors of M"^, where 
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k < d. Let K he a. cone which is generated by {^1,^2,...,^^}, in other words, 
K = {AiMi + X2U2 + . . . + XkUk, for some Aj > and i = 1,2,..., k}. Consider the 
parallelepiped S = jAiUi + A^tt ? + . . . + XkUk, < Aj < 1, z = 1, 2, . . . , k}. 



It is well-known 
K equals 



19971 ) that the generating function for the lattice points in 



k 



Thus, to derive a formula for arbitrary pointed cones one could decompose them 
into simple cones, via a trian gulation, and t hen apply the formula above and the 



inclusion-exclusion principle in 



2 



StanlevI ()l980l Proposition 1.2). Instead, Barvinok's 



idea is that it is more efficient to further decompose each simple cone into simple 
unimodular cones. A unimodular cone is a simple cone with generators {ui, . . . , Uk\ 
that form an integral basis for the lattice ]R{mi, . . . , Uk] H Z'^. Note that in this case 
the numerator of the formula has a single monomial, in other words, the parallelepiped 
has only one lattice point. 



4.1.1 Simple signed decompositions 

We now focus our attention on how the cone decomposition is done. To decompose 
a cone into simple cones the first step is to do a triangulation [triangulation of a 
cone K in dimension d is a collection of d- dimensional simple cones such that their 
union is K, their interiors are disjoint, and any pair of them intersect in a (possibly 
empty) common face). T here are efficient algorithrns, when the dimen sion is fixed, to 



carry a triangulation (see 



Aurenhammer and KleinI (|2000| ): lLeel (|l997|) for details). In 



LattE we use the well-known Delaunay triangulation which we compute via a convex 
hull calculation. The idea is to "lift" the rays of the cone into a higher dimensional 
paraboloid by adding a new coordinate which is the sum of the squares of the other 
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coordinates, take the lower convex hull of the lifted points, an d then ' 



project" back 



those simple facets. We use Fukuda's implementation in CDD fiFukudal . 120011 ) of this 
lift-and-project algorithm. This is not the only choice of triangulation, and definitely 
not the smallest one. In Section 14.41 we discuss some situations when the choice of 
triangulation in fact gives a better rational function. 
In principle, one could at this point list the points of the funda mental para l lelepip ed, 



for exampl e, using a fast 



iilbert bases code such as 4ti2 (Hemmecke 



NORMALIZ ( Bruns and Kock 



2002) or 



20011 ). and then use formula (*) for a general simple 
cone. Theoretically this is bad because the number of lattice points in the paral- 
lelepiped is exponentially large already for fixed dimension. In practice, this can 
often be done and in some situations is useful. Barvinok instead decomposes each 
simple cone as a (signed) sum of simple unimodular cones. To be more formal, for a 
set A C W^, the indicator function [A\ : M'^ ^ M of A is defined as 



[A]{x) 



1 if X e A, 

if a; ^ A. 

We want to express the indicator function of a simple cone as an integer linear combi- 
nation of the indicator functions of unimodular simple cones. There is a nice valuation 
from the algebra of indicator f unctio ns of polyhedra to the field of rational functions 



( Barvinok and Pommersheim 



19991), and many of its properties can be used in the 



calculation. For example, the valuation is zero when the polyhedron contains a line. 



Theorem 4.3. \Barvinok and Pommersheim . 



199yl. Theorem, 3.1) There is a valu- 



ation f from the algebra of indicator functions of rational polyhedra into the field of 
multivariate rational functions such that for any polyhedron P, f{[P]) = J^aePnz^ ^"^ ■ 

Therefore once we have a unimodular cone decomposition, the rational generating 
function of the original cone is a signed sum of "short" rational functions. Next we 
focus on how to decompose a simple cone into unimodular cones. 
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Let Ui,U2, ■ ■ ■ ,Ud be linearly independent integral vectors which generate a simple 
cone K. We denote the index of K by ind{K) which tells how far K is from being 
unimodular. That is, md{K) = \ det(ui|n2| . . . \ud)\ which is the volume of the par- 
allelepiped spanned hj ui,U2, . . . , u^. It is also equal to the number of lattice points 
inside the half-open parallelepiped. K is unimodular if and only if the index of K is 
1. Now we discuss how we implemented the following key result of Barvinok: 



Theorem 4.4. JiBarvinok and Pommersheim . 



199^. Theorem 4-2) Fix the dimension 



d. Then, there exists a polynomial time algorithm with a given rational polyhedral cone 
K C W^ , which computes unimodular cones Ki, i G / = {1,2, . . . ,/}, and numbers 
ej G { — 1,1} such that 

i&I 

Let K be a rational pointed simple cone. Consider the closed parallelepiped 

r = {aiUi + a2U2 + . . . + adUd : \aj\ < {ind{K))~d, j = 1, 2, . . . , d}. 

Note that this parallelepiped F is centrally sym metric and one c an show that the 
volume of F is 2'^. Minkowski's First Theorem (jSchrijveii . 119861 1 guarantees that 
because F C M'' is a centrally symmetric convex body with volume > 2"^, there exists 
a non-zero lattice point w inside of F. We will use w to build the decom position. 



Dver and Kannan 



We ne ed to find w explicitly. We take essentially the approach suggested by 

( 19931 ). We require a subroutine that computes the shortest vector in a lattice. For 

fixed dimension this can be don e in polynomial time using lattice basis reduction (this 



follows trivially from 



Schriivej (Corollary 6.4b ll986l page 72)). It i s worth obse rving 



that when the dimension is not fixed the problem becomes NP-hard (lAitaiL 



use the basis reducti on algorithm of Lenstra, Lenstra, and Lovasz (JGrotschel et al. . 



199.4 



Schriiver 



199fih. We 



1986^ to find a short vector. Given A, an integral d^ d matrix whose 



columns generate a lattice, LLL's algorithm outputs A' , a new d^d matrix, spanning 
the same lattice generated by A. The column vectors of A', m'^,M2, ...,«[;, are short 
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and nearly orthogonal to each other, and each u'^ is an approximatio n of the s 



lortest 



vector in the lattice, in terms of Euclidean length. It is well-known (jSchriiver 
that there exists a unique unimodular matrix U suc h that AU = A'. 



19m 



The method proposed in 



Dyer and KannanI ()l993| ) to find w is the following: Let 



A = (ui\u2\ ■ ■ ■ \ud), where the Ui are the rays of the simple cone we wish to decompose. 
Compute the reduced basis of A~^ using the LLL algorithm. Let A' = {u[\u2\ . . . \u'j) 
be the reduced basis of A^^. Dyer and Kannan observed that we can find the smallest 
vector with respect to the l°° norm by searching over all linear integral combinations of 
the column vectors of A' with small coefficients. We call this search the enumeration 
step. This enumeration step can be performed in polynomial time in fixed dimension. 
We will briefiy describe the enumeration step. First we introduce some notation. Let 
u[, . . . ,u'^he linearly independent integral vectors in Z''. Let || ■ II2 be the P norm 
and let || ■ ||oo be the infinity norm. 

We will need to recall the Gram-Schmidt process that computes a set of orthogonal 
vectors u*, 1 < j < d, from independent vectors u'j, 1 < j < d. In particular we need 
some values from this process. The vectors m*, and real numbers ^j^k-, 1 "^ k < j < d 
are computed from u'j by the recursive process: 






i-i 



«j = «j 



'j-'^f^j^kul, 2<j<d 



fe=i 



l^j,k 



u'j ■ K 



U 



*||2 
k\\2 



Letting Wi := u*/||n,*||2, there exists real numbers Ui{j) such that 



d 



Note that Ui{j) = fii.j\\u*\\2 for 1 < k < j < d and Ui{i) = \\u*\\2. These Ui{j) will be 
used below. Let L{u'i, . . . , u'^) be the lattice generated by m'^^, . . . , u'^. Then we denote 
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Lj{u[, . . . , u'^) be the projection of L{u[, . . . , u'^) orthogonal to the vector space Vj 

spanned hj u[, . . . , u'y 

Now we are ready to describe the process of the enumeration step. Let A be a shortest 

vector in the lattice spanned by A' with respect to the /°° norm. Then we can write 

A as an integral linear combination of columns of A' . Let A = X]i=i '^i'^'v where 

a = (q!i, . . . , ad) G Z"^. The goal is to find some finite set T <Z U^ such that a ^ T 

and the cardinality of T is polynomial size in fixed dimension. T will be contained 

inside a certain parallelepiped. Then we can search A by enumerating all lattice points 

inside T. 

Suppose A' = {u[\ . . . \u'^) form the reduced basis obtained by LLL algorithm. Let 

m := min{j : Uj{j) > Ui{l)} — 1. Now we will apply the inequalities 

lkl|oo< ||x||2, (4.2) 



\X\\2 



< v^||x||oo- (4.3) 



We are going to prove that a shortest vector of L{u[, . . . , u'^) is a shortest vector of 
L{u[, . . . , u'^) with respect to the /°° norm. Any vector in L{u[, . . . , u'^)\L{u[, . . . , u'^) 
must have /^ norm at least mi(1). Since mi(1) = ||wi||2 it must have /°° norm at least 
ni(l) which is at least the l°° norm of u'l by ()4.2|) . 

We will show how to construct T. Let y = X]i=i '^i'^'i be a candidate for a shortest 
vector with respect to the l°° norm . App lying the fact that |ajnj(z)| < \\y\\2 (using the 



same trick as on page 423. 



Kannan 



1 lloo 



( 19871 )). we have \aiUi(i)\/\/d < \\y\\2/Vd ^ il^ 
for any candidate vector for a shortest vector with respect to the l°° norm. Therefore, 
we have 



\ai\ui{{)/\Q < ||%||oo < ui{l). 

From this 

< vdui{l)/ui{i) for i = 1, . . . , m, 



Oii 
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which defines a parallelepiped in the variables a^ such that, 

Q:= {aeR'^ : -Vdui{l)/ui(i) < at < Vdui{l) / Ui{i) for i = 1, . . . , d}. 

Finally we set T := Q n Z'^. 

Now we are going to show that Q contains polynomially many lattice points. For 
each ttj, there exist at most 1 + 2\/dui{l) / Ui{i) candidates. So the total number of 
candidates is 

m 

J](l + 2v^ni(l)/n,(z)). 
With the fact that mi(1) > Ui{i) (by the definition of m), we have 

m m 

J](l + 2v^n(l)i/n,(z)) < S'^rf'^/^ J](ni(l)/n,(z)). 
We derive the following from Minkowski's theorem 

u,{ir<{2mr/'det{L{u[,...,u'J), 

m 

det(L(ui,...,M^)) = JjMi(i). 
Therefore, we have YY^=i{ui{\) / Ui{i)) < (2m)'"/^. This implies that 

m 

J](l + 2v^n(l)i/ni(z)) < (3c?)^ 

i=l 

which is a constant if we fix d. With this method, we can compute a shortest vector A 

with respect to the /°° norm in polynomial time in fixed dimension by the enumeration 

step. 

After we compute A in polynomial time in fixed dimension, we know that there exists 

a unique unimodular matrix U such that A' = A~^U . Minkowski's theorem for the 

l°° norm implies that for the non-singular matrix A' , there exists a non-zero integral 

vector ^: s uch t hat ||A||oo = ||^'-2||oo < | det(y4')|^/'^. See statement 23 in page 81 in 



Schrijveii ()l98fi[ ). We can set 
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< I det(A')|'/'' = I det{A~^U)\^/'' = I det{A~^) det([/)|^/'' 



= I det{A-')\'/'^ = I det{A)\-'^'^ = |ind(i^)r'/'='. 

Since A~^ and A' span the same lattice, there exists an integral vector w ^ Mf^ such 
that A = A~^w. Then, we have 

w = AX. 



'■) 



Note that ty is a non-zero integral vector which is a linear integer combination of the 
generators Ui of the cone K with possibly negative coefficients., and with coefficients 
at most |ind(i<')|~^/'^. Therefore, we have found a non-zero integral vector w E V. 
In LattE, we try to avoid the enumeration step because it is very costly. Instead, 
we choose A to be the shortest of the columns in A' . This may not be the smallest 
vector, but for practical purposes, it often decreases the index |ind(ii')| just like for a 
shortest vector. Experimentally we have observed that we rarely use the enumeration 
step. 
In the next step of the algorithm, for z = 1, 2, . . . , d, we set 

Ki = conejui, Us, . . . , n^.i, w, n^+i, . . . , n^}. 

Now, we have to show that for each i, ind{Ki) is smaller than ind(J'i'). Let w = 
J2i=i (^iUi- Then, we have 

ind(ii'i) = \ det{{ui\u2\ . . . \ui_i\w\ui+i\ . . . \ud))\ 

= |Q;i||det((Mi|M2| . . . \ui-i\ui\ui+i\ . . . \ud))\ 
= \ai\md{K) < (ind(A'))'^ . 
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There is one more technical condition that w needs to satisfy. This is that w and 



Ui, . . . ,Ud belong to an open half-space (jBarvinok . 



1994 Lemma 5.2). This is easy to 



achieve as either the w we found or —w satisfy this condition. We can now decompose 
the original cone K into cones Ki for i = 1,2, . . . ,d, oi smaller index, [K] = ^ ±[fCj]. 
This sum of indicator functions carries signs which depend on the position of w with 
respect to the interior or exterior of K. We iterate this process until Ki becomes 
a unimodular cone for i = 1,2,..., d. For i mplem enting Barvinok's decomposition 



of cones, we use the package NTL by 



Shoupl (|2003l ) to compute the reduced basis of 



a cone and to compute with matrices and determinants. All our calculations were 
done in exact long integer arithmetic using the routines integrated in NTL. Here is the 
pseudo-code of the algorithm and an example. 

Algorithm 4.5. (Barvinok's Decomposition of a Simple Cone) 
Input: A simple cone K = cone{ui,U2, ■ ■ ■ ,Ud\ given by its generators. 
Output: A list of unimodular cones and numbers ej as in Theorem\4.4 



Set two queues Uni and NonUni. 
if K is unimodular 

then Uni = Uni U {K}. 
else NonUni = NonUni U {K}. 
while Non Uni is not empty do 

Take a cone K G NonUni and set A = {ui, U2, ■ ■ ■ , Ud) 

to be a matrix whose columns are the rays of K. 
Compute the smallest vector A in the lattice, 

with respect to l°° , which is spanned by the column vectors of A~^. 
Find a non-zero integral vector z such that A = A~^z. 
if vectors z, U\, U2, ■ ■ ■ ,Ud are in an open half plane 
then set z := z. 
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else set z := —z. 
for % = 1, 2, . . . , d do 

set Ki = cone{ui, . . ., Mi_i, z, Ui+i, . . . ,Ud} 
and set Ai = (mi, . . . , Mj_i, z, Ui+i, ...,Ud). 
for z = 1, 2, . . . , (i do 

if det(Aj) and (\ei{A) have the same sign 

then assign e^^. = e^- 
else ek, = -ex- 
for i = 1,2, . . . ,d do 
if f^j is unimodular 

then [/nz = Uni yj{Ki}. 
else NonUni = NonUni U{Ki}. 
return all elements in Uni. 

It is very important to remark tliat, in principle, one also needs to keep track of 
lower dimensional cones present in the decomposition for the purpose of writing the 
inclusion-exclusion formula of the generating function f{K, z). For example in Figure 
14.21 we have counted a ray twice, and thus it needs to be removed. 




Figure 4.2: Contribution of lower dimensional cones 



But t his is actually not necessary thanks to a Brian's polarization trick (JBarvinok and Pommersheim 
19991 Remark 4.3): Let K* be the dual cone to K. Apply the iterative procedure 



above to K* instead of K, ignoring the lower dimensional cones. This can be done 
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because once we polarize the result back, the contribution of the lower dimensional 
cones is zero with respect to the valuation that assigns to an indicator funct i on its 



generating function counting the lattice points (jBarvinok and Pommersheim 



199S, 



Corollary 2.8). In the current implementation of LattE we do the following: 

1. Find the vertices of the polytope and their defining supporting cones. 

2. Compute the polar cone to each of the cones. 

3. Apply the Barvinok decomposition to each of the polars. 

4. Polarize back the cones to obtain a decomposition, into full- dimensional unimod- 
ular cones, of the original supporting cones. 

5. Recover the generating function of each cone and, by Brion's theorem, of the 
whole polytope. 

Here is an example of how we carry out the decomposition. 
Example 4.6. Let K be a cone generated by (2, 7)"^ and (1,0)"^. Let 



A 



2 1 
7 



Then, we have det{A) = —7 and 




The reduced basis A' of A ^ and the unimodular matrix U for the transformation from 



i f \ / 1 

A to A' are: A' = \ , and U = \ \ . By enumerating the column 

vectors, we can verify that (^, y)^ is the smallest vector with respect to l°° in the 
lattice generated by the column vectors of A~^. So, we have z = (1,0)"^. Then, we 
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have two cones: 

^ 2 \ / 1 

and I 

7 1 y \ ^ "^ 

The second cone is unimodular of index —1 which is the same sign as the determinant 
of A. Thus, Uni = Uni U { I }, and assign to it e = 1. The first cone has 

determinant 2. 5*0, we assign e = — 1. Since the first cone is not unimodular, we have 

.20, 
NonUm = Nonllni U{ | | }. Set 

7 1 



A 



2 

7 1 



Then, we have det{A) = 2 and 



A-i = I ^ M , A' = I ^ '^] andU=\ ^ ^ 

Since X = (|, ^)"^ is the smallest vector with respect to /°°, we have z = {1, 3)'^. So, 

we get two cones: 

^ 2 1 \ /^ 1 

and I 

7 3 y y 3 1 

The first matrix has negative determinant which is not the same sign as the determi- 
nant of its parent matrix A. Since eA = —1, we assign to the first cone e = 1 and the 
second one has positive determinant, so we assign to it e = 1. Since both of them are 
unimodular, we take them into Uni and since Nonllni is empty, we end while loop 
and print all elements in Uni. 
This gives a full decomposition: 
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Figure 4.3: Example of Barvinok's decomposition 



From the previous example, we notice that the determinant of each cone gets much 
smaller in each step. This is not an accident as Theorem 14.41 guarantees that the 
cardinality of the index set / of cones in the decomposition is bounded polynomially 
in terms of the determinant of the input matrix. We have looked experimentally 
at how many levels of iteration are necessary to carry out the decomposition. We 
observed experimentally that it often grows linearly with the dimension. We tested 
two kinds of instances. We used random square matrices whose entries are between 
and 9, thinking of their columns as the generators of a cone centered at the origin. 
We tested from 2x2 matrices all the way to 8 x 8 matrices, and we tested fifteen 
random square matrices for each dimension. We show the results in Table 14.11 For 
computation, we used a 1 GHz Pentium PC machine running Red Hat Linux. 



The seco nd set of examp 



matrices (jSchriiver 



es comes from the Birkhoff polytope Bn of doubly stochastic 



19861 ). Each vertex of the polytope is a permutation matrix which 
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Dimension 


Height of tree 


# of cones 


determinant 


Time (seconds) 


2 


1.33 


2.53 


11.53 





3 


2.87 


12.47 


55.73 


0.005 


4 


3.87 


65.67 


274.667 


0.153 


5 


5.87 


859.4 


3875.87 


0.25 


6 


7.47 


10308 


19310.4 


3.67 


7 


8.53 


91029.4 


72986.3 


41.61 


8 


10.67 


2482647.533 


1133094.733 


2554.478 



Table 4.1: Averages of 15 random matrices for computational experiences 



Dimension 


# of vertices 


# of unimodular cones at a vertex cone 


Time (seconds) 


S3 = 4 


6 


3 


0.05 


54 = 9 


24 


16 


0.15 


B5 = 16 


120 


125 


0.5 


Be = 25 


720 


1296 


7.8 



Table 4.2: The numbers of unimodular cones for the Birkhoff polytopes 



is a 0/1 matrix whose column sums and row sums are all 1 (jSchrijverl . |l2Sa)- We 
decompose the cone with vertex at the origin and whose rays are the n\ permutation 
matrices. The results are reported in Table W^ 



4.1.2 From cones to rational functions and counting 

Once we decompose all cones into simple unimodular cones, it is easy to find the gen- 
erating function attached to the ith cone Ki. In the denominator there is a product 
of binomials of the form (1 — z^'j) where Bij is the jth ray of the cone Ki. Thus the 
denominator is the polynomial 11(1 ~ z^'^). How about the numerator? The cone Ki 
is unimodular, thus it must have a single monomial z"^', corresponding to the unique 
lattice point inside the fundamental parallelepiped of Ki. Remember that the vertex 



CHAPTER 4. Experimental results: development of LattE 64 

of Ki is one of the vertices of our input polytope. If that vertex v has all integer coor- 
dinates then Ai = v, or else v can be written as a linear combination ^ ^jBij where 
all the Aj are rational numbers and can be found by solving a system of equations 
(remember the Bij form a vector space basis for M.'^). The unique lattice point inside 



the p arallelepiped of the cone Ki is simply ^ [Aj] Bij (JBarvinok and Pommersheim 



1999L Lemma 4.1). 

Brion's theorem says the sum of the rational functions coming from the unimodular 
cones at the vertices is a polynomial with one monomial per lattice point inside the 
input polytope. One might think that to compute the number of lattice points inside 
of a given convex polyhedron, one could directly substitute the value of 1 at each of 
the variables. Unfortunately, (1,1,..., 1) is a singularity of all the rational functions. 



Instead we discuss the method use d in LattE to compute this va 



from that presented by Barvinok (JBarvinok and Pommersheim 



ue, which is different 



19991). The typical 



generating function of lattice points inside a unimodular cone forms: 



EH] 



^A, 



n(i-2^'0' 

where z"' is monomial in d variables, each A^ (cone vertex) and B^j (a generator of 
cone i) are integer vectors of length d, i ranges over all cones given, j ranges over 
the generators of cone i, and E[i\ is 1 or -1. Adding these rational functions and 
simplifying would yield the polynomial function of the lattice point of the polytope. 
Now this is practically impossible as the number of monomials is too large. But 
calculating the number of monomials in this polynomial is equivalent to evaluating 
the limit as Zi goes to 1 for all i. We begin by finding an integer vector A and 
making the substitution Zi — ^ t^^ . This is with the intention of obtaining a univariate 
polynomial. To do this, A must be picked such that there is no zero denominator 
in any cone expression, i.e. no dot product of A with a Bij can be zero. Barvinok 
showed that such a A can be picked in polynomial time by choosing points on the 
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moment curve. Unfortunately, this method yields large values in the entries of A. 
Instead we try random vectors with small integer entries, allowing small increments 
if necessary, until we find A. Since we are essentially trying to avoid a measure zero 
set, this process terminates very quickly in practice. 

After substitution, we have expressions of the form ±t^^ / Yli^ ~ t^'^), where Ni and 
Dij are integers. Notice that this substitution followed by summing these expressions 
yields the same polynomial as would result from first summing and then substituting. 
This follows from the fact that we can take Laurent series expansions, and the sum 
of Laurent series is equal to the Laurent series of the sum of the original expressions. 
Also, note that we have the following identity: 

# of cones j^. 

After substitution we have the following univariate (Laurent) polynomial such that: 

# of cones .y, 

aePnZ'* i=i ^^^ ' 

With the purpose of avoiding large exponents in the numerators, we factor out a 
power of t, say f^. Now we need to evaluate the sum of these expressions at t = 1, 
but we cannot evaluate these expressions directly at t = 1 because each has a pole 
there. Consider the Laurent expansion of the sum of these expressions about t = 1. 
The expansion must evaluate at t = 1 to the finite number X^^gpn^d 1. It is a Taylor 
expansion and its value at t = 1 is simply the constant coefficient. If we expand each 
expression about t = 1 individually and add them up, it will yield the same result as 
adding the expressions and then expanding (again the sum of Laurent expansions is 
the Laurent expansion of the sum of the expressions). Thus, to obtain the constant 
coefficient of the sum, we add up the constant coefficients of the expansions about 
t = 1 of each summand. Computationally, this is accomplished by substituting 
t = s + 1 and expanding about s = via a polynomial division. Summing up 
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the constant coefficients with proper accounting for E[i] and proper decimal accuracy 
yields the desired result: the number of lattice points in the polytope. Before the 
substitution t = s + 1 we rewrite each rational function in the sum (recall f^ was 
factored to keep exponents small); 



t^'-'^ ^ „,. , t^^ 



E^[^]n7T37l^ = E^'[ 



involves in such a way that /},•• > for all i,j. This requires that the powers of t at 
each numerator to be modified, and the sign E[i] is also adjusted to E'[i]. Then the 
substitution t = s + 1 yields 

^ n((i + s)°--i) 

where it is evident that, in each summand, the pole s = has an order equal to the 
number of factors in the denominator. This is the same as the number of rays in the 
corresponding cone and we denote this number by d. 

Thus the summand for cone i can be rewritten as E'[i]s~'^Pi{s)/Qi{s) where Pi{s) = 
(1 + s)^' and Qi{s) = Y[ {{^ + s)^'^ — 1) / s) . Pi{s)/Qi{s) is a Taylor polynomial whose 
s'^ coefficient is the contribution we are looking for (after accounting for the sign E'[i] 
of course). The coefficients of the quotient Pi{s)/Qi{s) can be obtained recursively 
as follows: Let Qi{s) = bo + biS + b2s'^ + . . . and -Pj(s) = ag + ais + a2S^ + . . . and let 
p4^ = Co + cis + C2S^ + . . .. Therefore, we want to obtain Cd which is the coefficient 
of the constant term of Pi/Qi- So, how do we obtain q from Qi{s) and Pi{s)7 We 
obtain this by the following recurrence relation: 



_ ao 

bo 

Ck = 7-(«fc - biCk-i - b2Ck-2 - ... - 6fcCo) for A; = 1, 2, ... . 
In order to obtain q, only the coefficients Oq, ai, . . . , a^ and b^, bi, . . . ,bd are required. 
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Example 4.7. (A triangle). Let us consider three points in 2 dimensions such that 
Vi = (0, 1), V2 = (1,0), and V3 = (0,0). Then, the convex hull ofVi, V2, and V3 is a 
triangle in 2 dimensions. We want to compute the number of lattice points by using 
the residue theorem. Let Ki be the vertex cone at Vi for i = 1,2,3. Then, we have 
the rational functions: 



^^'^^' ^"' y^^ = (1 - y-4l - xy~^y ^^"^^ ^"' '^^ = (l-x-^)(l-x^^y) ' 

f{K3,{x,y)) ^ 



{l-x){l-yy 
We choose a vector A such that the inner products of A and the generators of Ki are 

not equal to zero. We choose A = (1, —1) in this example. Then, reduce multivariate 

to univariate with X, so that we have: 



t-' 



^^^^' ^^ - {i-t){i-t^y ^^^^' ^^ - (i-t-i)(i-t-2)' /(^^3, t) - ^Y^TxT^FT)- 

We want to have all the denominators to have positive exponents. We simplify them 
in order to eliminate negative exponents in the denominators with simple algebra. 
Then, we have: 



t^ 



f{K,,t) = ——-——-, /(K2,t) = ,,, J{K3.t) 



(l-t)(l-t2)'^^ ^' ^ (l_t)(l_t2)'^V ^'"^ (l_t)(l_t)' 

We factor out t^^ from each rational function, so that we obtain: 



fiK,,t) = — — ^-— -,/(i^2,t) = — — — — -,/(i^3,t) 



:i-t)(l-t2)'^^ "' ^ (l_t)(l_t2)'^V -'"^ (l_t)(l-t)' 



We substitute t = s + 1 and simplify them to the form 






,,,, , 1 ,,r^ , 1 + 5S + lOs^ + 10s3 + 5^4 + ^5 

s^[2 + s) s'^[2 + s) 
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f{K3,s) 



,2 



Now we use the recurrence relation to obtain the coefficient of the constant terms. 
Then, for f{Ki,s), we have C2 = |. For f{K2,s), we have C2 = ^. For f{K^,s), 
we have C2 = — 1. Thus, if we sum up all these coefficients, we have 3, which is the 
number of lattice points in this triangle, q 

LattE produces the sum of rational functions which converges to the generating func- 
tion of the lattice points of an input polytope. This generating function is a multivari- 
ate polynomial of finite degree. As we saw in Subsection 14. 1 . 21 it is possible to count 
the number of lattice points without expanding the rational functions into the sum of 
monomials. Suppose that instead of wanting to know the number of lattice points we 



simply wish to decide whether there is one lattice point inside the po 



The i nteger feasibility problem is an important and difficult problem (jAardal et al. 



1998; 



Schriiver 



ytope or not 



19861 ). Obviously, one can simply compute the residues and then if 
the number of lattice points is non-zero, clearly, the polytope has lattice points. 
Before we end our description of LattE, we must comment on how we deal with 
polytopes that are not full-dimensional (e.g. transportation polytopes). Given the 
lower-dimensional polytope P = {x G M"' : Ax = a, Bx < h} with the d x n matrix 
A of full row-rank, we will use the equations to transform P into a polytope Q = 
{x G M""'* : Cx < c} in fewer variables, whose integer points are in one-to-one 
correspondence to the integer points of P. This second polytope will be the input to 
the main part of LattE. The main idea of this transformation is to find the general 
integer solution x = a;o + X]r=i -^«5'« ^° ^^ ~ ^ ^^^ ^^ substitute it into the inequalities 
Bx < b, giving a new system Cx < c in n — d variables Ai, . . . , Xn~d- 



It is known that the general i nteger so 



ution Ax = a can be found via the Hermite 



normal form H = {R\0) of A (jSchriiveii . 119861 ). Here, i? is a lower-triangular matrix 



and H = AU for some unimodular matrix U . Moreover, as A is supposed to have full 
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row-rank, i? is a non-singular d x d matrix. Let f/i be the matrix consist of the first 
d columns of U and U2 consisting of the remaining n — d columns of U. Now we have 
AUi = R and AU2 = and the columns of U2 give the generators {gi, . . . ,gn-d} of 
the integer null-space of A. Thus, it remains to determine a special integer solution 
xq to Ax = a. 

To do this, first find an integer solution y^ to Hy = {R\0)y = a, which is easy due 
to the triangular structure of R. With xq = Ui/o, we get Axq = AUyo = Hy^ = a 
and have found all pieces of the general integer solution x = Xq + X]r=i '^idi ^° 
{x G Z" : Ax = a}. 



4.2 Computational experience and performance for counting 



One can download LattE via a web page |www. math ■ucdavis.edu/~latte| You can 
also find there the files of all the experiments presented in this section. At the moment 
we have been able to handle polytopes of dimension 30 and several thousands vertices. 
It is known that the theoretical upper bound of the number of unimodular cones is 
2dh -^j^gj^g Ji = I °s og ■ ~ °^ °^ — I and where D is the volume of the fundamental 



1994^ . If we fix the dimension this upper 



parallelepiped of the input cone (jBarvinok , 
bound becomes polynomial time. Unfortunately, if we do not fix the dimension, this 
upper bound becomes exponential. In practice this might be costly and some families 
of polytopes have large numbers of unimodular cones. The cross polytope family, 
for instance, has many unimodular cones and behaves b adly. For exa mple, for the 



20011), LattE took 



cross polytope in 6 dimensions, with crossG.ine input file (jFukudal 

147.63 seconds to finish computing. The number of lattice points of this polytope is 

obviously 13. A lso, for the cross polytope in 8 dimensions, with crossS.ine input file 



(Pukuda 



200ll ). LattE took 85311.3 seconds to finish computing, even though this 



polytope has only 16 vertices and the number of lattice points of this polytope is 17. 
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For all computations, we used a 1 GHz Pentium PC machine running Red Hat Linux. 

Here is a short description of how to use LattE. 

For computations involving a polytope P described by a system of inequalities Ax < b, 



"tjj, 



and b G Z™, the LattE readable input file would be as 



where A e Z""""^, A -- 
follows: 

m d+1 
b -A 

EXAMPLE. Let P = {{x,y) -.x <l,y <l,x + y <l,x>0,y>0}. Thus 

/ 1 0\ /l\ 

1 1 

A= 11 , 6= 1 

-10 

V -1/ Vo/ 

and the LattE input file would be as such: 

5 3 

1-10 
10-1 
1 -1 -1 
1 
1 

In LattE, polytopes are represented by linear constraints, i.e. equalities or inequal- 
ities. By default a constraint is an inequality of type ax < b unless we specify, by 
using a single additional line, the line numbers of constraints that are linear equalities. 
EXAMPLE. Let P be as in the previous example, but require x + y = 1 instead of 
X + y < 1, thus, P = {{x,y) : x < l,y < l,x + y = l,x > 0,y >0}. Then the LattE 
input file that describes P would be as such: 
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5 3 
1-10 
10-1 
1 -1 -1 
1 

1 
linearity 1 3 

The last line states that among the 5 inequalities one is to be considered an equality, 
the third one. 

For bigger examples it quickly becomes cumbersome to state all nonnegativity con- 
straints for the variables one by one. Instead, you may use another short-hand. 
EXAMPLE. Let P be as in the previous example, then the LattE input file that 
describes P could also be described as such: 

3 3 

1-10 

10-1 

1 -1 -1 
linearity 1 3 
nonnegative 2 12 

The last line states that there are two nonnegativity constraints and that the first 
and second variables are required to be nonnegative. NOTE that the first line reads 
"3 3" and not "5 3" as above! 

We now report on computations with convex rational polytopes. We used a 1 GHz 
Pentium PC machine running Red Hat Linux. We begin with the class of multiway 
contingency tables. A d-table of size (rii, . . . , n^) is an array of non-negative integers 
V = ivii,...,ij, 1 < "^j < nj. For < m < d, an m-marginal of v is any of the 
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(^) possible ?7i-tables obtained by summing the entries over all but m indices. For 
instance, if (vij^k) is a 3-table then its 0-marginal is f+,+,+ = Xir=i Xl?ii J2^Li'^i,j,k, 
its 1-marginals are (wj^+^+) = (X]?=i J2^Li'^'i,j,k) and likewise (t'+j,+), (f+,+,fc), and its 
2-marginals are (t'ij,+) = (Z]fcli^ij,fc) and likewise (wi,+,fc), (^'+,j,fc)- 
Such tables appear naturally in statistics and operations research under various names 
such as multi-way contingency tables, or tabular data. We consider the table counting 
problem: given a prescribed collection of marginals, how many d-tables are there that 
share these marginals? Table counting has several applications in statistical analysis, 
i n particular for independence testing, and has been the focus of much research (see 



( Diaconis and Gangolli . 



19951 ) and the extensive list of references therein). Given a 
specified collection of marginals for d-tables of size (rii, . . . , Ud) (possibly together with 
specified lower and upper bounds on some of the table entries) the associated multi- 
index transportation polytope is the set of all non-negative real valued arrays satisfying 
these marginals and entry bounds. The counting problem can be formulated as that 
of counting the number of integer points in the associated multi-index transportation 
polytope. We begin with a small example of a three-dimensional table of format 
2x3x3 given below. The data displa yed in Table 14.3 have been extracted from 



the 1990 decennial census and is used in (JFienberg et al. 



20011 ). For the 2-marginals 



implied by these data we get the answer of 441 in less than a second. 

We present now an example of a 3 x 3 x 3 table with fairly large 2-marginals. They are 

displayed in Table U3] LattE took only 19.67 seconds of CPU time. The number of 

lattice points inside of this polytope is 2249847900174017152559270967589010977293. 

Next we present an example of a 3 x 3 x 4 table with large 2-marginals. The 2- 

marginals are displayed in Table 1^31 The CPU time for this example was 44 minutes 

42.22 seconds. The number of lattice points inside of this polytope is 

4091700129572445106288079361219676736812805058988286839062994. 

The next family of examples are some hard knapsack-type problems. Suppose we have 
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Gender = Male 

Income Level 



Race 


< $10,000 


> $10000 and < $25000 


> $25000 


Total 


White 


96 


72 


161 


329 


Black 


10 


7 


6 


23 


Chinese 


1 


1 


2 


4 


Total 


107 


80 


169 


356 



Gender = Female 

Income Level 



Race 


< $10,000 


> $10000 and < $25000 


> $25000 


Total 


White 


186 


127 


51 


364 


Black 


11 


7 


3 


21 


Chinese 





1 





1 


Total 


197 


135 


54 


386 



Table 4.3: Three-way cross-classification of gender, race, and income for a selected U.S. census tract. 
Source: 1990 Census Public Use Microdata Files. 

a set of positive relatively prime integers {ai,a2, . . • ,0,^}. Denote by a the vector 
(oi, 02, . . . , Crf). Consider the following problem: does there exist a non-negative 
integral vector x satisfying ax = an for some positive integer qq? We take several 



examples from (jAardal et al. 



2002aJ) which have been found to be extremely hard to 



solve by commercial quality branch-and-bound software. This is very surprising since 
the number of variables is at most 10. It is not very difficult to see that if the right- 
hand-side value Oo is large enough, the equation will surely have a non-negative integer 
solution. The Frobenius number for a knapsack problem is the largest value On such 



that the knapsack problem is infeasible. Aardal and Lenstr a (|Aardal et al. 



solved them using the reformulation in (jAardal et al 



2002a|) 



1998|). Their method works 
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164424 


324745 


127239 


262784 


601074 


9369116 


149654 


7618489 


1736281 




163445 


49395 


403568 


1151824 


767866 


8313284 


1609500 


6331023 


1563901 




184032 


123585 


269245 


886393 


6722333 


935582 


1854344 


302366 


9075926 



Table 4.4: 2-Marginals for the 3 x 3 x 3 example. 



273510 


273510 


273510 


191457 


273510 


273510 


547020 


191457 


273510 


547020 


273510 


191457 



464967 


273510 


273510 


547020 


273510 


464967 


410265 


601722 


273510 



273510 


273510 


273510 


410265 


547020 


136755 


547020 


136755 


410265 


191457 


191457 


191457 



Table 4.5: 2-Marginals for the 3x3x4 example. 
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significantly better than branch-and-bound using CPLEX 6.5. Here we demonstrate 



that our implementation of Barvinok's algorithm is fairly fa st and, on t 



le order 



of seconds, we resolved the first 15 problems in Table 1 of (jAardal et al.L |2002a|) 
and verified all are infeasible except probQ, where there is a mistake. The vector 
(3480, 1, 4, 4, 1, 0, 0, 0, 0, 0) is a solution to the right-hand-side 13385099. In fact, using 
LattE we know that the exact number of solution s is 838908602000. Fo r comparison 



we named the problems exactly as in Table 1 of (jAardal et al.L |2002aJ). We present 
our results in Table 14.61 It is very interesting to know the number of lattice points 
if we add 1 to the Frobenius number for each problem. In Table 14. 7^ we find the 
number of solutions if we add 1 to the Frobenius number on each of the (infeasible) 
problems. The speed is practically the same as in the previous case. In fact the speed 
is the same regardless of the right-hand-side value Qq. 

Aheady counting the lattice points of large width convex polygons is a non-trivial task 
if one uses brute-force enumeration (e.g. list one by one the points in a bounding box 
of the polygon and see whether it is inside the polygon). Here we experiment with very 
large convex almost regular n-gons. Regular n-gons cannot have rational coordinates, 
but we can approximate them to any desired accuracy by rational polygons. In the 
following experiment we take regular n-gons, from n = 5 to n = 12 centered at the 
origin (these have only a handful of lattice points). We take a truncation of the 
coordinates up to 3, 9, and 15 decimal digits, then we multiply by a large enough 
power of 10 to make those vertex coordinates integral and we count the number of 
lattice points in the dilation. All experiments take less than a second. 
The next two sets of examples are families that have been studied quite extensively 
in the literature and provide us with a test for speed. In the first case we deal with 
two-way contingency tables. The polytope defined by a two-way contingency table is 
called the transportation polytope. We present the results in Table 14.21 The second 
family consists of flow polytopes for the complete 4-vertex and the complete 5- vertex 







Frobenius # 


Time (m, s) 


cuwwl 


12223 


12224 


36674 


61119 


85569 








89643481 


0.55s 


cuww2 


12228 


36679 


36682 


48908 


61139 


73365 






89716838 


1.78s 


cuww3 


12137 


24269 


36405 


36407 


48545 


60683 






58925134 


1.27s 


cuww4 


13211 


13212 


39638 


52844 


66060 


79268 


92482 




104723595 


2.042s 


cuww5 


13429 


26850 


26855 


40280 


40281 


53711 


53714 


67141 


45094583 


16.05s 


prol 


25067 


49300 


49717 


62124 


87608 


88025 


113673 


119169 


33367335 


47.07s 


prob2 


11948 


23330 


30635 


44197 


92754 


123389 


136951 


140745 


14215206 


lm0.58s 


prob3 


39559 


61679 


79625 


99658 


133404 


137071 


159757 


173977 


58424799 


lm28.3s 


prob4 


48709 


55893 


62177 


65919 


86271 


87692 


102881 


109765 


60575665 


59.04s 


prob5 


28637 


48198 


80330 


91980 


102221 


135518 


165564 


176049 


62442884 


lm41.78s 


prob6 


20601 


40429 


40429 


45415 


53725 


61919 


64470 


69340 78539 95043 


22382774 


3m45.86s 


prob7 


18902 


26720 


34538 


34868 


49201 


49531 


65167 


66800 84069 137179 


27267751 


2m57.64s 


probS 


17035 


45529 


48317 


48506 


86120 


100178 


112464 


115819 125128 129688 


21733990 


8m29.78s 


problO 


45276 


70778 


86911 


92634 


97839 


125941 


134269 


141033 147279 153525 


106925261 


4m24.67s 



o 

> 

H 



X 



< 
a> 

o' 

xs 
B 



M 



Table 4.6: Infeasible knapsack problems. 



as 
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problem 


RHS 


# of lattice points. 


cuwwl 


89643482 


1 


cuww2 


89716839 


1 


cuww3 


58925135 


2 


cuww4 


104723596 


1 


cuww5 


45094584 


1 


probl 


33367336 


859202692 


prob2 


14215207 


2047107 


prob3 


58424800 


35534465752 


pro4 


60575666 


63192351 


pro5 


62442885 


21789552314 


pro6 


22382775 


218842 


pro7 


27267752 


4198350819898 


pro8 


21733991 


6743959 


pro 10 


106925262 


102401413506276371 



Table 4.7: The number of lattice points if we add 1 to the Frobenius number. 





10^ (seconds) 


10'-* (seconds) 


10^^ (seconds) 


5gon 


2371673(0.136) 


2377641287748905186(0.191) 


2377641290737895844565559026875(0.289) 


6gon 


2596011(0.153) 


2598076216000000011(0.193) 


2598076211353321000000000000081(0.267) 


7gon 


2737110(0.175) 


2736410188781217941(0.318) 


2736410188638105174143840143912(0.584s) 


8gon 


2820021(0.202) 


2828427120000000081(0.331) 


2828427124746200000000000000201(0.761) 


9gon 


2892811(0.212) 


2892544245156317460(0.461) 


2892544243589428566861745742966(0.813) 


lOgon 


2931453(0.221) 


2938926257659276211(0.380) 


2938926261462380264188126524437(0.702) 


llgon 


2974213(0.236) 


2973524496796366173(0.745) 


2973524496005786351949189500315(1.858) 


12gon 


2997201(0.255) 


3000000004942878881(0.466) 


3000000000000005419798779796241(0.696) 



Table 4.8: The numbers of the approximated regular polygons. We show the number of lattice 
points in different dilation factors (powers of ten) and time of computation. 
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tournaments (directed complete graphs). Consider the directed complete graph Ki 
for / G N and / > 3. We assign a number to each node of the graph. Then, we 
orient the arcs from the node of smaller index to the node of bigger index. Let A^ 
be the node set of the complete graph Ki, let Wi be a weight assigned to node i for 
i = 1, 2, . . . , /, and let A be the arc set of Ki. Then, we have the following constraints, 
with as many variables as arcs: 



Wi, 



(j,i)arc entersj {jj)arc has taih 

Xij >0 V(i,j) G A 

These equalities and inequalities define a polytope and this polytope is the special 
case of a flow polytope. The results for the complete graphs 7^4 and K^, with different 
weight vectors, are shown in Tables E. 101 and l4. Ill respectivelv. 
hese two families of poly t opes have been studied by several authors 



f Baldoni-Silva et al. 



2003; 



Beck 



2003; 



De Loera and Sturmfels . 



2001; 



Mount 



20001 ) 



and thus are good for testing the performance of LattE. We used several examples 



of transportation polytopes, as presented in t he table below. In genera. 



at comparable performance to the software of (JBaldoni-Silva et al. . 



2003; 



LattE runs 



Beck 



20031) 



for generic vectors (a, b) but is slower for degenerate inputs (those that do not give 

a simple polytope). The reason seems to be that at each non-simplex vertex LattE 

needs to triangulate each cone which takes considerable time in problems of high 

dimension. 

The following experiment is from JDiaconis and Sturmfeld (|l998h . This is from an 



actual data in German survey. For 2,262 German citizens they asked the following 
questions: if you order the following items, how you order them? 

1. Maintain order. 



2. Give people more say. 



Table 4.9: Testing for 4x4 transportation polytopes. 



O 
> 

H 



Margins 


# of lattice points 


Time (seconds) 


[220, 215, 93, 64], 
[108, 286, 71, 127] 


1225914276768514 


1.048 


[109, 127, 69, 109], 
[119, 86, 108, 101] 


993810896945891 


1.785 


[72, 67, 47, 96], 
[70, 70, 51, 91] 


25387360604030 


1.648 


[179909, 258827, 224919, 61909], 
[190019, 90636, 276208, 168701] 


13571026063401838164668296635065899923152079 


1.954 


[229623, 259723, 132135, 310952], 
[279858, 170568, 297181, 184826] 


646911395459296645200004000804003243371154862 


1.765 


[249961, 232006, 150459, 200438], 
[222515, 130701, 278288, 201360] 


319720249600111437887229255487847845310463475 


1.854 


[140648, 296472, 130724, 309173], 
[240223, 223149, 218763, 194882] 


322773560821008856417270275950599107061263625 


1.903 


[65205, 189726, 233525, 170004], 
[137007, 87762, 274082, 159609] 


6977523720740024241056075121611021139576010 


1.541 


[251746, 282451, 184389, 194442], 
[146933, 239421, 267665, 259009] 


861316343280649049593236132155039190682027614 


1.880 


[138498, 166344, 187928, 186942], 
[228834, 138788, 189477, 122613] 


63313101414342827754566531364533378588986467 


1.973 


[20812723, 17301709, 21133745, 27679151], 
[28343568, 18410455, 19751834, 20421471] 


665711555567792389878908993624629379187969880179721169068827951 


2.917 


[15663004, 19519372, 14722354, 22325971]. 
[17617837, 25267522, 20146447, 9198895] 


63292704423041655080293971395348848807454253204720526472462015 


3.161 


[13070380, 18156451, 13365203, 20567424], 
[12268303, 20733257, 17743591, 14414307] 


43075357146173570402117291685601604830544643769252831337342557 


2.990 



X 



(TI 



< 
a> 

o' 

xs 
B 



M 



CO 
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Weights on nodes 


# of lattice points 


Time (seconds) 


[-6, -8, 5, 9] 


223 


0.288 


[-9, -11, 12, 8] 


330 


0.286 


[-1000, -1, 1000, 1] 


3002 


0.287 


[-4383, 886, 2777, 720] 


785528058 


0.287 


[-4907, -2218, 3812, 3313] 


20673947895 


0.288 


[-2569, -3820, 1108, 5281] 


14100406254 


0.282 


[-3842, -3945, 6744, 1043] 


1906669380 


0.281 


[-47896, -30744, 46242, 32398] 


19470466783680 


0.282 


[-54915, -97874, 64165, 88624] 


106036300535520 


0.281 


[-69295, -62008, 28678, 102625] 


179777378508547 


0.282 


[-3125352, -6257694, 926385, 8456661] 


34441480172695101274 


0.509 


[-2738090, -6701290, 190120, 9249260] 


28493245103068590026 


0.463 


[-6860556, -1727289, 934435, 7653410] 


91608082255943644656 


0.503 



Table 4.10: Testing for the complete graph K4 



Weights on nodes 


# of lattice points 


sees 


[-12, -8, 9, 7, 4] 


14805 


0.319 


[-125, -50, 75, 33, 67] 


6950747024 


0.325 


[-763, -41, 227, 89, 488] 


222850218035543 


0.325 


[-11675, -88765, 25610, 64072, 10758] 


563408416219655157542748 


0.319 


[-78301, -24083, 22274, 19326, 60784] 


1 108629405 144880240444547243 


0.336 


[-52541, -88985, 1112, 55665, 84749] 


3997121684242603301444265332 


0.331 


[-71799, -80011, 86060, 39543, 26207] 


160949617742851302259767600 


0.316 


[-45617, -46855, 24133, 54922, 13417] 


15711217216898158096466094 


0.285 


[-54915, -97874, 64165, 86807, 1817] 


102815492358112722152328 


0.277 


[-69295, -62008, 28678, 88725, 13900] 


65348330279808617817420057 


0.288 


[-8959393, -2901013, 85873, 533630, 11240903] 


6817997013081449330251623043931489475270 


0.555 


[-2738090, -6701290, 190120, 347397, 8901863] 


277145720781272784955528774814729345461 


0.599 


[-6860556, -1727289, 934435, 818368, 6835042] 


710305971948234346520365668331191134724 


0.478 



Table 4.11: Testing for the complete graph K5. Time is given in seconds 
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1234 


137 


2134 


48 


3124 


330 


4123 


21 


1243 


29 


2143 


23 


3142 


294 


4132 


30 


1324 


309 


2314 


61 


3214 


117 


4213 


29 


1342 


255 


2341 


55 


3241 


69 


4231 


52 


1423 


52 


2413 


33 


3412 


70 


4312 


35 


1432 


93 


2431 


39 


3421 


34 


4321 


27 


Total 


875 




279 




914 




194 



Table 4.12: Permutation 5*4 problem from 



Diaconis and Sturmfela ( 1998 ) 



875 


279 


914 


194 


746 


433 


742 


341 


345 


773 


419 


725 


296 


777 


187 


1002 



Table 4.13: Marginal conditions for the permutation problem from 



Diaconis and Sturmfelj (|l99S ) 



3. Fight rising prices. 



4. Protect freedom of speech. 

Then we have the data in Table 14.121 

We take sum for the number of people who picked j G {1, 2, 3, 4} for the ith order for 

all j = 1, 2, 3, 4 and for all i = 1,2, 3, 4. Then we have the following condition given 

in Table mH 

Then we want to compute how many functions / such that f : S4 ^ N and 

E.es, /(^) = 2262. 

The solution to this problem is: Total number of functions is 
11606690287805167142987310121 and CPU Time is 523.12 sec. 
The last experiment in this section is from Ian Dinwoodie. Ian Dinwoodie commu- 
nicated to us the problem of counting all 7 x 7 contingency tables whose entries are 
nonnegative integers Xi, with diagonal entries multiplied by a constant as presented 



CHAPTER 4. 



Experimental results: development of LattE 



82 



2x1 


X2 


2:3 


X4 


X5 


X6 


X7 


205 


X2 


2X8 


Xg 


a;io 


Xll 


a;i2 


a;i3 


600 


X3 


Xg 


2X14 


a;i5 


2^16 


Xi7 


a:i8 


61 


X4 


Xio 


Xl5 


2X19 


X20 


a;2i 


3:22 


17 


X5 


Xll 


a;i6 


X20 


2X23 


a;24 


X2b 


11 


Xq 


2^12 


Xi7 


a;2i 


2^24 


2X26 


X2I 


152 


Xl 


Xl3. 


xia 


2;22 


X2b 


a;27 


2X28 


36 


205 


600 


61 


17 


11 


152 


36 


1082 



Table 4.14: The conditions for retinoblastoma RBl-VNTR genotype data from the Ceph database. 



in Table 14. 141 The row sums and column sums of the entries are given there too. 
Using LattE we obtained the exact answer 8813835312287964978894. 



4.3 New Ehrhart (quasi-)polynoniials 

Given a rational polytope P C W^, the function 

Zp(t):=#(tPnZ^), 



19771 ) and has 



for a positive integer t, was first studied by E. Ehrhart (jEhrhart 

received a lot of attention in combinatorics. It is known to be a polynomial when all 

vertices of P are integral and it is a quasi-polynomial for arbitrary rati onal polytopes . 



1997 . 



It is called the Ehrhart quasi-polynomial in honor of its discoverer (jStanlev 
Chapter 4). A function / : N — *> C is a quasi-polynomial if there is an integer A^ > 
and polynomials fo, . . . , /n-i such that f{s) = fi{s) if s = z mod A^. The integer 
A^ is called a quasi-period of /. Therefore, by counting the number of lattice points 
for sufficiently many dilations of a rational polytope, we can interpolate its Ehrhart 
quasi-polynomial. 

Using LattE, Maple, and interpolation, we have calculated the Ehrhart polynomials 
and quasi-polynomials for polytopes that are slices or nice truncations of the unit 
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d-cuhe. To the best of our knowledge these values were not known before. For 
example, the 24-cell polytope centered at the origin with smallest integer coordinates 
has Ehrhart polynomial i24_cea{s) = 8s'^ + ^ + 8s'^ + ^ + 1. In Table KT^ we 
see the Ehrhart polynomials for the hypersimplices A{n,k). They are defined as 
the slice of the n-cube by the hyperplane of equation ^ Xj = /c with k < n. Note 
that A{n,k) = A{n,n — k) because of the symmetries of the regular cube. The 
hypersimplices form one of the most famous f amilies of 0/1-po l ytope s. It is known 



20011 ). This means 



that hypersimplices are compressed polytopes (jOhsugi and Hibi . 

that their Ehrhart polynomials can be recovered from the /-vectors of any of their 

reverse lexicographic triangulations. Instead, we recovered them explicitly for the 

first time using LattE and interpolation. 

We also have the Ehrhart quasi-polynomials of some truncated unit cubes. 

Proposition 4.8. The Ehrhart quasi-polynomial for the truncated unit cube in Figure 
\4-4i where its vertices are at 1/3 and 2/3 of the way along edges of the cube, is given 
by: 



'•tru^cubeiy^ ) 



77^3 I 235^ I 19 
81 "•" 9 "•" 9 



+ 1 if s = mod 3, 
^ + ^-|f-f z/s^l mod 3, 
- |i ifs = 2 mod 3. 



TTsf _|_ 6l£ 
81 ~'~ 27 

77s3 



81 ~'~ 27 



29s 

27 




Figure 4.4: The truncated cube. 
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Proposition 4.9. The Ehrhart quasi-polynomial for the cuboctahedron (Figure \J~^ 



is: 



'^tru_cube2\S ) 



^ + 2s^ + f + lifs = mod 2, 
^ + ^-f-fz/s = l mod 2. 




Figure 4.5: The cuboctahedron. 

Proposition 4.10. The Ehrhart quasi-polynomial for the truncated regular simplex, 
where the vertices are at 1/3 and 2/3 of the way along the simplex edges (see Figure 
\4-(>\ l, is given by: 



''tru_simplex\S j 



2'is-^ 



9 



13s 
9 



1 if s = mod 3, 



1f + ¥ + t + S^/^-2 mod3. 



(1, 0, 1) 




Figure 4.6: The truncated simplex. 
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4.4 Computational results via Homogenized Barvinok's al- 
gorithm 

We have demonstrated the practical relevance of Barvinok's cone decomposition ap- 
proach for counting lattice points and deriving formulas. Several other algorithms are 
available to carry out the same kind of enumeration. It is important to implement 
them all in the same computer system for comparison of performance and to corrob- 
orate that the answers are correct. Some problems are solvable by some methods but 
not by others. 
We would lik e to remind a rea der that one of the bottle necks of the original Barvinok 



algorithm in (JBarvinok 



1994 ) is the fact that a polytope may have too many vertices. 
Originally, we visit all vertices of the input polytope to compute Barvinok's short 
rational function and this can be costly in terms of computational time. For example, 
the well-known polytope of semi-magic cubes in the 4x4x4 case has over two million 
vertices, but only 64 linear inequalities describe the polytope. Algorithm 12. 101 shows 
that Homogenized Barvinok's algorithm works with only a single cone. In this section, 
we will show some practical results with Homogenized Barvinok's algorithm. 
A normal semigroup S is the intersection of the lattice Z*^ with a rational convex 
polyhedral cone in M.'^. Each pointed affine semigroup S <Z U^ can be graded. This 
means that there is a linear map deg : 5 ^ N with deg{x) = if and only if a; = 0. 
Given a pointed graded affine semigroup, we define Sr to be the set of elements with 
degree r, i.e. Sr = {x & S : deg{x) = r}. The Hilbert series of S is the formal power 
series Hs(t) = J^'kLo i^i^r)^^, where #(5*^) is the cardinality of Sr- Algebraically, this 
is just the Hilbert series of the semigroup ring C[S]. It is a well-known property that 
Hs is represented by a rational function of the form 

Q{t) 



[I - t'^){i - 1'^) . . . {1 - 1 



Sd' 
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wher e Q(t) is a polynomial of degree less than Si + ■ ■ ■ + s^ (see Chapter 4 ([Stanley 



19971)). The first challenge is to compute the Hilbert series of magic cubes. Sev- 



eral other met 



fl Ahmed et ah. 



lods had been tried to compute the Hilbert series explicitly (see 



20031 ) for references). One of the most well-known challenges was 
that of counting the 5x5 magic squares of magic sum n. Similarly several authors 
had tried before to compute the Hilbert series of the 3x3x3x3 magic cubes. It 
is not difficult to see this is equivalent to determining an Ehrhart series. Using Algo- 
rithmic^! we finally present the solution, which had been inaccessible using Grobner 
bases methods. For comparison, the reader familiar with Grobner bases computations 
should be aware that the 5x5 magic squares problem required a computation of a 
Grobner bases of a toric ideal of a matrix A with 25 rows and over 4828 columns. Our 
attempts to handle this problem with CoCoA and Macaulay2 were unsuccessful. We 
now give the numerator and then the denominator of the rational functions computed 
with the software LattE: 

Theorem 4.11. 

The generating function X]n>o/('^)^" Z^'" ^^^ number f{n) of 5 x 5 magic squares of 
magic sum n is given by the rational function p{t) / q{t) with numerator 

p{t) = t^6 + 28 1^5 + 639 f^ + 11050 f^ + 136266 1^^ _|_ 1255833 f^ + 9120009 t™ + 54389347 1'^^ + 
274778754 1'^^ + 1204206107 1'^^ + 4663304831 1^^ + 16193751710 1®^ + 51030919095 1®"* + 
147368813970 1'^^ + 393197605792 t^^ + 975980866856 t'^^ + 2266977091533 t^° + 4952467350549 4^9 + 
10220353765317 t^s + 20000425620982 4^7 + 37238997469701 1^*^ + 66164771134709 t^s + 112476891429452 t^* + 
183365550921732 t^^ + 287269293973236 t^^ + 433289919534912 t^i + 630230390692834 t^o 
+ 885291593024017 t*^ + 1202550133880678 t^^ + 1581424159799051 i"'' + 2015395674628040 t""^ + 
2491275358809867 t*^ + 2989255690350053 t"*"* + 3483898479782320 t"*^ + 3946056312532923 t*^ + 
4345559454316341 1"! + 4654344257066635 t^° + 4849590327731195 i^^ + 4916398325176454 t^^ + 
4849590327731195 t^'^ + 4654344257066635 t^'^ + 4345559454316341 1^^ + 3946056312532923 t^** + 
3483898479782320 t^^ + 2989255690350053 t^^ + 2491275358809867 t^^ + 2015395674628040 t^° + 
1581424159799051 i^^ + 1202550133880678 i^** + 885291593024017 t'^'' + 630230390692834 t^e + 433289919534912 t^s + 
287269293973236 t^* + 183365550921732 t^^ + 112476891429452 t^a + 66164771134709 t^i + 37238997469701 1'^° + 
20000425620982 t^^ + 10220353765317 1^** + 4952467350549 t^'^ + 2266977091533 t^*^ + 975980866856 t^^ + 
393197605792*14 + 147368813970 1^^ + 51030919095 i^^ + 16193751710 i" +4663304831*10 + 1204206107*9 + 
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274778754 1^ + 54389347 f + 9120009 1« + 1255833 1^ + 136266 f* + 11050 1^ + 639 1^ + 28 1 + 1 and denominator 

q(t) = (t2 - 1) 1° {t^+t+iy {f ~ 1)' (t6 +t^ +1) {t'- +f^ +t^ +t+l)Ul~ tf {t^ + 1)". 

The generating function Yln>ofi^)^"^ /^'^ ^^^ number f{n) o/ 3x3x3x3 magic 
cubes with magic sum n is given the rational function r{t)/s{t) where 

r(t) = t54 + 150 {51 + 5837^1** + 63127 1''^ + 331124 1"'^ + 1056374 t^s + 2326380 1^*^ + 3842273 i^^ + 5055138 4^" + 
5512456 127 + 5055138*2** + 3842273*^1 + 2326380**8 + 1056374**5 + 331124**2 ^ 63127*^ + 5837*"^ + 150 *3 + 1 
and 

S{t) = (*3 + 1)* (**2 + t9 + t6 + *3 + 1) (1 - t3)9 (46 + *3 + 1) . 

4.5 Computational results via the BBS algorithm and the 
digging algorithm 



In this se ction we report our expe rience solving hard knapsack problems from 



Aardal et al 



()2002a^ : ICornueiols et all (|l997|). See Table WTEi for the data used here. Their form 
is maximizec ■ x subject to ax = b,x >0,x E Z'*, where b E 1j and where a 'Elf' with 
gcd(ai, . . . , Orf) = 1. For the cost vector c, we choose the first d components of the vec- 
tor (213,-1928,-11111,-2345,9123,-12834,-123,122331,0,0). We compared the 
original digging algorithm, the single cone digging algorithm, and the BBS algorithm, 
which are implemented in LattE (available at |www . math . ucdavis . edu/~latte]) , with 
CPLEX version 6.6. The computations were done on a 1 GHz Pentium PC running 
Red Hat Linux. Table 14.171 provides the optimal values and an optimal solution for 
each problem. As it turns out, there is exactly one optimal solution for each problem. 
With one exception, CPLEX 6.6. could not solve the given problems. Note that 
whenever the digging algorithm found the optimal value, it did so much faster than 
the BBS algorithm. This is interesting, as the worst-case complexity for the digging 
algorithm is exponential even for fixed dimension, while the BBS has polynomial 
complexity in fixed dimension. The digging algorithm fails to find a solution for 
problems cuww2, cuww3, and cuww5. What happens is that the expansion step 
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Problem 


a 


b 


cuwwl 


12223 


12224 


36674 


61119 


85569 












89643482 


cuww2 


1222S 


36679 


36682 


48908 


61139 


73365 










89716839 


cuww3 


12137 


24269 


36405 


36407 


48545 


60683 










58925135 


cuww4 


13211 


13212 


39638 


52844 


66060 


79268 


D2482 








104723596 


cuwwS 


13429 


26850 


26855 


40280 


40281 


53711 


53714 


67141 






45094584 


probl 


25067 


49300 


49717 


62124 


87608 


88025 


113673 


119169 






33367336 


pi:ob2 


11948 


23330 


30635 


44197 


92754 


123389 


136951 


140745 






14215207 


probS 


39559 


61679 


79625 


99658 


133404 


137071 


159757 


173977 






58424800 


prob4 


48709 


55893 


62177 


65919 


86271 


87692 


102881 


109765 






60575666 


probS 


28637 


48198 


80330 


91980 


102221 


135518 


165564 


176049 






62442885 


probG 


20601 


40429 


40429 


45415 


53725 


61919 


64470 


69340 


78539 


95043 


22382775 


probT 


18902 


26720 


34538 


34868 


49201 


49531 


65167 


66800 


84069 


137179 


27267752 


probS 


17035 


45529 


48317 


48506 


86120 


100178 


112464 


115819 


12512S 


129688 


21733991 


probQ 


3719 


20289 


29067 


60517 


64354 


65633 


76D6D 


102024 


106036 


119930 


13385100 


problO 


45276 


70778 


86911 


92634 


97839 


125941 


134260 


141033 


147279 


153525 


106925262 



Table 4.16: knapsack problems. 

becomes costly when more coefficients have to be computed. In these three examples, 
we computed coefficients for more than 2,500,000, 400,000, and 100,000 powers of t; 
all turning out to be 0. The Digging algorithm is slower than CPLEX in problem 
prob9 because duri ng the execution of Barvinok's unimodular cone decomposition (see 



pages 15 and 16 of 



Barvinok and PommersheimI (|1999|)) more than 160,000 cones are 



generated, leading to an enormous rational function for f{P;t). Moreover, for prob9 
more than 3,500 coefficients turned out to be 0, before a non-zero leading coefficient 
was detected. Finally, in problems cuwwl, cuww3, prob2, probS, prob4, prob6, and 
probS, no digging was necessary at all, that is, Lasserre's condition did not fail here. 
For all other problems, Lasserre's condition did fail and digging steps were necessary 
to find the ffist non- vanishing coefficient in the expansion of /(P; t). 



Problem 


Value 


Solution 


Runtime for 
Digging (Original) 


Runtime for 
Digging (S. Cone) 


Runtime for 
BBS 


Runtime for 
CPLEX 6.6 


cuwwl 


1562142 


[7334 0] 


0.4 sec. 


0.17 sec. 


414 sec. 


> 1.5h (OM) 


cuww2 


-4713321 


[3 2445 0] 


> 3.51i 


> 3.5h 


6,600 sec. 


> 0.75b (OM) 


cuww3 


1034115 


[4855 0] 


1.4 sec. 


0.24 sec. 


6,126 sec. 


> 0.75b (OM) 


cuww4 


-29355262 


[0 2642 0] 


> 1.51i 


> 1.5h 


38,511 sec. 


> 0.75b (OM) 


cuww5 


-3246082 


[1 1678 10 0] 


> 1.51i 


147.63 sec. 


> sob 


> 0.75b (OM) 


probl 


9257735 


[966 5 1 74] 


51.4 sec. 


18.55 sec. 


> 3h 


> Ih (OM) 


prob2 


3471390 


[853 2 4 27] 


24.8 sec. 


6.07 sec 


> lOh 


> 0.75b (OM) 


prob3 


21291722 


[708 2 1 173] 


48.2 sec. 


9.03 sec. 


> 12h 


> 1.5h (OM) 


prob4 


6765166 


[1113 7 54] 


34.2 sec. 


9.61 sec. 


> 5h 


> 1.5h (OM) 


prob5 


12903963 


[1540 12 103] 


34.5 sec. 


9.94 sec. 


> 5h 


> 1.5h (OM) 


prob6 


2645069 


[1012 1 1 1 20 0] 


143.2 sec. 


19.21 sec. 


> 4h 


> 2h (OM) 


prob7 


22915859 


[782 1 1 186 0] 


142.3 sec. 


12.84 sec. 


> 4h 


> Ih (OM) 


probS 


3546296 


[1 385 1 1 35 0] 


469.9 sec. 


49.21 sec. 


> 3.5h 


> 2.5h (OM) 


prob9 


15507976 


[31 11 1 1 127 0] 


1,408.2 sec. 


283.34 sec. 


> lib 


4.7 sec. 


problO 


47946931 


[0 705 1 1 403 0] 


250.6 sec. 


29.28 sec. 


> lib 


> Ih (OM) 



o 

> 

H 



X 



(TI 



< 
a> 

o' 

xs 
B 



M 



Table 4.17: Optimal values, optimal solutions, and running times for eacb problem. 0M:= Out of memory. 
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problem 


Original 
Digging (A) 


Original 
Digging (B) 


Single Cone 
Digging (A) 


Single Cone 
Digging (B) 


cuww 1 


110 





25 





cuww 2 


386 


> 2,500,000 


79 


> 2,500,000 


cuww 3 


346 





49 





cuww 4 


364 


> 400,000 


51 


> 400,000 


cuww 5 


2,514 


> 100,000 


453 


578,535 


prob 1 


10,618 


74,150 


1,665 


74,150 


prob 2 


6,244 





806 





prob 3 


12,972 





2,151 





prob 4 


9,732 





1,367 





prob 5 


8,414 


1 


2,336 


1 


prob 6 


26,448 


5 


3,418 


5 


prob 7 


20,192 





2,015 





prob 8 


62,044 





6,523 





prob 9 


162,035 


3,558 


45,017 


3,510 


prob 10 


38,638 


256 


5,128 


256 



Table 4.18: Data for the digging algorithm. A := munber of unimodular cones and B := number of 
digging levels 
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problem 


(C) 


(D) 


(E) 


euww 1 


125,562 


4,829.3 


26 


cuww 2 


2,216,554 


88,662.16 


25 


euww 3 


2,007,512 


80,300.48 


25 


cuww 4 


10,055,730 


402,229.2 


25 


cuww 5 


NA 


NA 


< 26 



Table 4.19: Data for the BBS algorithm. C := Total num. of unimodular cones, D := Average num. 
of unimodular cones per an iteration, E := Total number of iterations, and NA := not available. 



problem 


The optimal value 


An optimal solution 


(N) 


cuww 1 


1562142 


[7334 0] 




cuww 2 


-4713321 


[3 2445 0] 




cuww 3 


1034115 


[4855 0] 




cuww 4 


-29355262 


[0 2642 0] 




cuww 5 


-3246082 


[1 1678 10 0] 




prob 1 


9257735 


[966 5 1 74] 




prob 2 


3471390 


[853 2 4 27] 




prob 3 


21291722 


[708 2 1 173] 




prob 4 


6765166 


[1113 7 54] 




prob 5 


12903963 


[1540 12 103] 




prob 6 


2645069 


[1012 1 1 1 20 0] 




prob 7 


22915859 


[782 1 1 186 0] 




prob 8 


3546296 


[1 385 1 1 35 0] 




prob 9 


15507976 


[31 11 1 1 127 0] 




prob 10 


47946931 


[0 705 1 1 403 0] 





Table 4.20: The optimal value and an optimal solution for each problem. N := num. of solutions 
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Appendix A 



User manual of LattE 



A.l Introduction 



A. 1.1 What is LattE? 



The name "LattE" is an abbreviation for "Lattice point Enumeration." So what 
exactly does LattE do? The software's main function is to count the lattice points 
contained in convex polyhedra defined by linear equations and inequalities with in- 
teger coefficients. The polyhedra can be of any (reasonably small) dimension, and 
LattE use s an algorithm that runs in po l ynorn ial time for fixed dimension: Barvinok's 



algorithm 



Barvinok and PommersheimI (jl999f ). To learn more about the exact details 



of our i mplementation a r id algo rithmic techniques involved, the interested reader can 



consult 



De Loera et al. 



(2(md 



M) 



and the references listed therein. Here we give 
a rather short description of the mathematical objects used by LattE, Barvinok's 
Rational Functions: 

Given a convex polyhedron P = {u & Mf^ : Au < 6}, where A and h are integral, the 



fundamental object that we compute is a short representation of the infinite power 
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series: 



1 ™a2 ™"d 

Xo . . . x^ . 



f{P;x) — 2_^ ^\'-^2 

Here each lattice point is given by one monomial. Note that this can be a rather 
long sum, in fact for a polyhedral cone it can be infinite, but the good news is that 
it admits short representations. 

Example: Let P be the quadrangle with vertices Vi = (0, 0), V2 = (5, 0), V3 = (4, 2), 
and 1/4 = (0,2). 







































V4 










X 


















\ 


















\ 


V, 






X 





































/(P; X, y) = x^ + x*y + a;" + x'^y^ + yx^ + a;''^ + x^y'^ + yx"^ + x^ + x'^y' ^ + xy + x + xy'^ + y + I + u 



'he fundamental theorem of Barvinok (circa 1993, see 



Barvinok and Pommersheim 



says that you can write /(P; x) as a sum of short rational functions, in 
polynomial time when the dimension of the polyhedron is fixed. In our running 
example we easily see that the 16 monomial polynomial can be written as shorter 
rational function sum: 

f{P;x,y) = f{Kv^;x,y) + f{Kv^;x,y) + f{Kvi,x,y) + f{Kv^;x,y) 
where 
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fiKv,;x,y) = ^^^^^^ f{Kv,;x,y) = ^^_,%'[t$.-^) 
f[Kv3]x,y) = jYZ^:^Tyjz^:p2^ f\Kv^:,x,y) = (i„j,-i)(i_^) 

/(P; 1,1) = 16 

Counting the lattice points in convex polyhedra is a powerful tool which allows many 

applications in areas such as Combinatorics, Statistics, Optimization, and Number 

Theory. 

A. 1.2 What can LattE compute? 

In the following we list the operations that LattE vl.l can perform on bounded convex 
polyhedra (more commonly referred to as polytopes). For the reader's convenience, 
we already include the basic commands to actually do the tasks. Let us assume that 
a description of a polytope P is given in the file "fileName" (see Section IA.3I for 
format) and that a cost vector is specified in the file "fileName. cost" (needed for the 
optimization part, see Section [A. 31 for format). 
Tasks performed by LattE vl.l: 

1. Count the number of lattice points in P. 

./count fileName 

2. Count the number of lattice points in nP, the dilation of P by the integer factor 
n. 

./count dil n fileName 

3. Calculate a rational function that encodes the Ehrhart series associated with 
the polytope. By definition, the n-th coefficient in the Ehrhart series equals the 
number of lattice points in nP . For more det ails on Ehrhart counting functions 



see, for example. Chapter 4 of 



StanlevI (|l997h . 
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. / ehrhart f i 1 eName 

4. Calculate the first n + 1 terms of the Ehrhart series associated with the polytope. 

./ehrhart n fileName 

5. Maximize or minimize a given linear function of the lattice points in P. 

./maximize fileName 
./minimize fileName 

In addition to these basic functions, there are more specific calls to LattE. For example 
to use the homogenized Barvinok algorithm instead of the original one in order to 
count the lattice points. These details will be explained in Section I A. 41 
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A. 2 Downloading and Installing LattE 

LattE is downloadable from the following website: 

littp://www. math. ucdavis.edu/~latte/downloads/ 
Step 1: Create directory for LattE 

mkdir latte 

Step 2: Download "latte_vl.l.tar.gz" to directory "latte" 

Download ' 'latte_vl . 1 .tar .gz' ' from 

http://www.math.ucdavis.edu/~latte/downloads/ 
(If you have never downloaded a file from the internet: A click with your right mouse 
button onto the file name on the webpage should do the trick. In any case, if you do 
not succeed, ask your system administrator, a friend, or send us an email.) 
Step 3: Change to directory for "latte" 

cd latte 

Step 4: Unzip and untar the archive 

gunzip latte_vl . 1 .tar .gz 
tar xvf latte_vl . 1 .tar 

Step 5: Make "install" executable 

chmod 700 install 
Step 6: Install LattE 
./install 
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A. 3 Input Files 

A. 3.1 LattE Input Files 
Inequality Description 

For computations involving a polytope P described by a system of inequalities Ax < b, 
where A G Z"*^'^, A = {aij), and b G Z™, the LattE readable input file would be as 
follows: 

m d+1 
b -A 

EXAMPLE. Let P = {{x,y):x <l,y <l,x + y <l,x>0,y>0}. Thus 

/ 1 o\ /l\ 

1 1 

A= 11 , 6= 1 

-10 

V -1/ yo/ 

and the LattE input file would be as such: 



5 


3 




1 


-1 





1 





-1 


1 


-1 


-1 





1 











1 


Equat 


ons 



In LattE, polytopes are represented by linear constraints, i.e. equalities or inequal- 
ities. By default a constraint is an inequality of type ax < b unless we specify, by 
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using a single additional line, the line numbers of constraints that are linear equalities. 
EXAMPLE. Let P be as in the previous example, but require x + y = 1 instead of 
X + y < 1, thus, P = {{x,y) : x < l,y < l,x + y = l,x > 0,y > 0}. Then the LattE 
input file that describes P would be as such: 

5 3 
1-10 
10-1 
1 -1 -1 
1 

1 
linearity 1 3 

The last line states that among the 5 inequalities one is to be considered an equality, 
the third one. 

Nonnegativity Constraints 

For bigger examples it quickly becomes cumbersome to state all nonnegativity con- 
straints for the variables one by one. Instead, you may use another short-hand. 
EXAMPLE. Let P be as in the previous example, then the LattE input file that 
describes P could also be described as such: 

3 3 

1-10 

10-1 

1 -1 -1 
linearity 1 3 
nonnegative 2 12 

The last line states that there are two nonnegativity constraints and that the first 



APPENDIX. User manual of LattE 100 

and second variables are required to be nonnegative. NOTE that the first hne reads 
"3 3" and not "5 3" as above! 

Cost Vector 

The functions maximize and minimize solve the integer linear programs 

max{c'^x : X e P n Z"'} 

and 

mm{c^x : X e P n Z'^}. 

Besides a description of the polyhedron P, these functions need a linear objective 
function given by a certain cost vector c. If the polyhedron is given in the file "file- 
Name" 

4 4 

1-10 

10-10 

10 0-1 

1 -1 -1 -1 

linearity 1 4 

nonnegative 3 12 3 

the cost vector must be given in the file "fileName.cost", as for example in the fol- 
lowing three-dimensional problem: 

1 3 

2 4 7 

The first two entries state the size of a 1 x n matrix (encoding the cost vector), 
followed by the 1 x n matrix itself. Assuming that we call maximize, this whole data 
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encodes the integer program 

max{2xi + 4x2 + 7x3 : Xi + X2 + X3 = 1, Xi, X2, X3 G {0, 1}}. 

A.3.2 cdd Input Files 

In addition to the formats described above, LattE can also accept input files in 
standard cdd format. (See Subsection IA.4. ll for details on how to run LattE on a cdd 
input file.) Below is an example of cdd input that is readable into LattE. 

H-representation 

begin 

4 4 integer 

2-2 4-1 

3-2-2 3 

6 2-4-3 

12 2 1 

end 

It is important to note that LattE can only read integer input. Clearly, cdd's rational 
data files can be converted into integer files by multiplying by the right constants. In 
the packaged release of LattE we include a binary version of cdd. 



APPENDIX. User manual of LattE 102 

A. 4 Running LattE 

A. 4.1 Command Syntax 

The basic syntax to invoke the various functions of LattE is: 

./count fileName 
. / ehrhart f i 1 eName 
./maximize fileName 
./minimize fileName 

Note that the last two functions require a cost vector specified in the file "file- 
Name. cost" ! 

Additionally, a variety of options can be used. All options should be space-delimited 
in the command. 

One option that can be set in addition to the options given below is "cdd" which tells 
LattE to read its input from a cdd input file. Thus, the above invocations for cdd 
input files would be 

./count cdd fileName 
./ehrhart cdd fileName 
./maximize cdd fileName 
./minimize cdd fileName 

A. 4. 2 Counting 

• Count the number of lattice points in polytope P, where P is given in "fileName" . 

./count fileName 

• Count the number of lattice points in nP, the dilation of P by the integer factor 
n. 
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./count dil n fileName 

• Count the number of lattice points in the interior of the polytope P, where P is 
given in "fileName". 

./count int fileName 



• Use the homogenized Barvinok algorithm iDe Loera et al.l ()2003b[ ) to count the 
number of lattice points in the polytope P, where P is given in "fileName". Use 
if number of vertices of P is big compared to the number of constraints. 

./count homog fileName 

A. 4. 3 Ehrhart Series 

• Compute the Ehrhart series encoded as a rational function for the polytope given 
in "fileName". Writes the unsimplified rational function to file "fileName. rat". 

. / ehrhart f i 1 eName 

• Compute the Ehrhart series encoded as a rational function for the polytope given 
in "fileName" . NEEDS Maple for simplification of terms. Writes the simplified 
rational function to file "fileName. rat". 

./ehrhart simplify fileName 

• Compute the Taylor series expansion of Ehrhart generating function up to degree 
n for the polytope given in "fileName" . 

./ehrhart n fileName 
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A. 4. 4 Optimizing 

This functions NEEDS a cost vector specified in "fileNamexost"!!! 

• Maximizes/Minimizes given linear cost functio n over the lattice point s in the 



De Loera et al. 



(2003^)isused. 



polytope given in "fileName" . Digging algorithm 
Optimal point and optimal value is returned. 

./maximize fileName 
./minimize fileName 

Maximizes/Minimizes given linear cost function over the lattice points in the 
polytope given in "fileName". Binary search algorithm is used. Only optimal 
value is returned. 

./maximize bbs fileName 
./minimize bbs fileName 
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A. 5 A Brief Tutorial 

In this section we invite the reader to follow along a few examples that show how to 
use LattE and also how to counter-check results. 



A. 5.1 Counting Magic Squares 

Our first example deals with counting magic 4x4 squares. We call a 4 x 4 array 
of nonnegative numbers a magic square if the sums of the 4 entries along each row, 
along each column and along the two main diagonals equals the same number s, the 
magic constant. Let us start with counting magic 4x4 squares that have the magic 
constant 1. Associating variables Xi,...,Xiq with the 16 entries, the conditions of 
a magic 4x4 square of magic sum 1 can be encoded into the following input file 
"EXAMPLES/magic4x4" for LattE. 

10 17 

1-1-1-1-1000000000000 

10000-1-1-1-100000000 

100000000-1-1-1-10000 

1000000000000-1-1-1-1 

1-1000-1000-1000-1000 

10-1000-1000-1000-100 

100-1000-1000-1000-10 

1000-1000-1000-1000-1 

1-10000-10000-10000-1 

1000-100-100-100-1000 

linearity 10 123456789 10 

nonnegative 16 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

Now we simply invoke the counting function of LattE by typing: 



APPENDIX. User manual of LattE 106 

./count EXAMPLES/magic4x4 

The last couple of lines that LattE prints to the screen look as follows: 

Total Unimodular Cones: 418 

Maximum number of simplicial cones in memory at once: 27 

***** Total number of lattice points: 8 **** 

Computation done. 
Time: 1.24219 sec 

Therefore, there are exactly 8 magic 4x4 squares that have the magic constant 1. 
This is not yet impressive, as we could have done that by hand. Therefore, let us 
try and find the corresponding number for the magic constant 12. Since this problem 
is a dilation (by factor 12) of the original problem, we do not have to create a new 
file. Instead, we use the option "dil" to indicate that we want to count the number 
of lattice points of a dilation of the given polytope: 

./count dil 12 EXAMPLES/magic4x4 

The last couple of lines that LattE prints to the screen look as follows: 

Total Unimodular Cones: 418 

Maximum number of simplicial cones in memory at once: 27 

***** Total number of lattice points: 225351 **** 

Computation done. 
Time: 1.22656 sec 

Therefore, there are exactly 225351 magic 4x4 squares that have the magic constant 
12. (We would NOT want to do THAT one by hand, would we?!) 
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Here is some amazing observation: the running time of LattE is roughly the same for 
counting magic squares of sum 1 and of sum 12. This phenomenon is due to the fact 
that the main part of the computation, the creation of the generating function that 
encodes all lattice points in the polytope, is nearly identical in both cases. 
Although we may be already happy with these simple counting results, let us be a 
bit more ambitious and and let us find a counting formula that, for given magic sum 
s, returns the number of magic 4x4 squares that have the magic constant s. 
For this, simply type (note that LattE invokes Maple to simplify intermediate expres- 
sions): 

./ehrhart simplify EXAMPLES/magic4x4 

The last couple of lines that LattE prints to the screen looks as follows: 

Rational function written to EXAMPLES/magic4x4.rat 

Computation done. 
Time: 0.724609 sec 

We are informed that this call created a file "EXAMPLES/magic4x4.rat" containing 
the Ehrhart series as a rational function: 

(t~8+4*t"7+18*t~6+36*t~5+50*t"4+36*t~3+18*t~2+4*t+l)/(-l+t)~4/(-l+t~2)~4 

Now we could use Maple (or your favorite computer algebra software) to find a series 
expansion of this expression. 

ts + 4 * t^ + 18 * t6 + 36 * t5 ^ 50 * t^ + 36 * t^ + 18 * t^ + 4 ^ t + 1 

(-l+t)4(-l+t2)4 

= 1 + 8t^ + 48^2 + 200t^ + 675t^ + 1904t^ + 4736t^ + 10608t^ + 21925t^ + 
42328t^ + 77328t^° + 134680t^^ + 225351^^^ + 364000t^^ + 570368t^^ + 
869856t^^ + 0(t^^) 
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The summands 8t and 225351t^^ reconfirm our previous counts. 

Although this rational function encodes the full Ehrhart series, it is not always as 

easy to compute as for magic 4x4 squares. As it turns out, adding and simplifying 

rational functions, although in just one variable t, can be extremely costly due to the 

high powers in t and due to long integer coefficients that appear. 

However, even if we cannot compute the full Ehrhart series, we can at least try and 

find the first couple of terms of it. 

./ehrhart 15 EXAMPLES/magic4x4 

The last couple of lines that LattE prints to the screen look as follows: 

Memory Save Mode: Taylor Expansion: 

1 

8t"l 

48t~2 

200t"3 

675t"4 

1904t~5 

4736t~6 

10608t"7 

21925t"8 

42328t"9 

77328t"10 

134680t~ll 

225351t~12 

364000t~13 

570368t~14 

869856t~15 
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Computation done. 
Time: 1.83789 sec 

Again, our previous counts are reconfirmed. 

Nice, but the more terms we want to compute the more time-consuming this task 
becomes. Clearly, if we could find sufficiently many terms, we could compute the full 
Ehrhart series expansion in terms of a rational function by interpolation. 

A. 5. 2 Counting Lattice Points in the 24-Cell 

Our next example deals with a well-known combinatorial object, the 24-cell. Its 
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10 10 
2-1 1 1-1 
10 0-1 
2-1-1 1 -1 
10-100 
2-1-1 1 1 
2 -1 -1 -1 1 
2 -1 -1 -1 -1 
1-10 

Now we invoke the counting function of LattE by typing: 

. / c ount EXAMPLES / 24_ c ell 

The last couple of lines that LattE prints to the screen look as follows: 

Total Unimodular Cones: 240 

Maximum number of simplicial cones in memory at once: 30 

***** Total number of lattice points: 33 **** 

Computation done. 
Time: 0.429686 sec 

Therefore, there are exactly 33 lattice points in the 24-cell. We get the same result 
by using the homogenized Barvinok algorithm: 

./count homog EXAMPLES/24_cell 

The last couple of lines that LattE prints to the screen look as follows: 

Memory Save Mode: Taylor Expansion: 
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**** Total number of lattice points is: 33 **** 

Computation done. 
Time: 0.957031 sec 

But how many of these 33 points he in the interior of the 24-ceU? 

./count int EXAMPLES/24_cell 

The last couple of lines that LattE prints to the screen look as follows: 

Reading .ext file... 



***** Total number of lattice points: 1 **** 

Therefore, there only one of the 33 lattice points in the 24-cell lies in the interior. 

A. 5.3 Maximizing Over a Knapsack Polytope 



Finally , let us solve the problem "cuwwl" 



Cornueiols et al 



(JMd) 
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( 2nn3ht) . Its description is given in the file "EXAMPLES/cuwwl": 



1 6 

89643482 -12223 -12224 -36674 -61119 -85569 

linearity 1 1 

nonnegative 5 12 3 4 5 

The cost function can be found in the file "EXAMPLES/cuwwl. cost": 

1 5 

213 -1928 -11111 -2345 9123 

Now let us maximize this cost function over the given knapsack polytope . Note that 

(2OO3J) is used. 



by default, the digging algorithm as described in 



De Loera et al. 
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./maximize EXAMPLES /cuwwl 
The last couple of lines that LattE prints to the screen look as follows: 

Finished computing a rational function. 
Time: 0.158203 sec. 

There is one optimal solution. 

No digging. 

An optimal solution for [213 -1928 -11111 -2345 9123] is: [7334 0]. 

The projected down opt value is: 191928257104 

The optimal value is: 1562142. 

The gap is: 7995261.806 

Computation done. 

Time: 0.203124 sec. 

The solution (7334, 0, 0, 0, 0) is quickly found. Now let us try to find the optimal 
value again by a different algorithm, the binary search algorithm. 

./maximize bbs EXAMPLES/cuwwl 

The last couple of lines that LattE prints to the screen look as follows: 

Total of Iterations: 26 

The total number of unimodular cones : 125562 

The optimal value: 1562142 

The number of optimal solutions: 1 
Time: 0.042968 

Note that we get the same optimal value, but no optimal solution is provided. 
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